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ANALYSIS OF NUMERICAL INTEGRATION TECHNIQUES FOR 
REAL-TIME DIGITAL FLIGHT SIMULATION 
By John W. Wilson and George G. Steinmetz 
Langley Research Center 

SUMMARY 

Low-order numerical integration techniques are analyzed and established as ade- 
quate for digital simulation of man-in-the-loop nonaerodynamic rigid-body problems. In 
addition to low-order integration techniques, a Green's function approach is presented for 
the solution of the gyroscope equations. This approach yields an exact solution to the 
gyroscope output for the simulated forcing function with virtually no limitation on the 
integration interval size. Effects of roundoff error are shown to be marginal for the 
translational dynamics with solution rates on the order of 32 solutions per second on a 
floating-point machine with a 27 -bit fractional part. 

Two example problems were solved by using the techniques established. A Gemini- 
Agena eleven-degree-of-freedom problem was programed in FORTRAN and then modified 
to a Gemini-Agena elastically coupled system. Primary concern was given to minimizing 
the time required of the arithmetic unit. The results of these problems indicate the util- 
ity of the IBM 7094 class of computers for this type of problem since only a small frac- 
tion of the memory and arithmetic processing time was used (10 to 15 percent each). 

This study strongly supports the desirability of multiprograming to alleviate inefficiency 
on a more sophisticated machine. 


INTRODUCTION 

In the past few years, general-purpose digital computers have played an important 
role in meeting resolution and accuracy requirements of real-time simulation problems. 
(See ref. 1, p. R34.) However, the extensive use of digital computing equipment has been 
prohibited (ref. 1, pp. R29-R30) except in special applications. (See ref. 2.) In the 
early years of the 1960 decade, general-purpose digital computers of sufficient speed and 
storage were available at a cost that was competitive (at least for the larger simulation 
problems; ref. 1, p. R31, fig. 3) with general-purpose analog equipment. Even so, the 
general-purpose analog facility, which is a composite of several basic computers, could 
be tied together in different combinations to perform one large simulation or several 



smaller simulations; this arrangement provides an efficiency advantage that the general- 
purpose digital computer did not yet enjoy. The capacity of the digital computer was 
determined by the largest problem to be solved and was totally dedicated to the solution 
of a smaller problem; thus, a large amount of idle central processor time and core 
remained unused. This fact led many simulation laboratories to develop software capa- 
ble of doing more than one simulation on a "one computer complex either simultaneously 
or alternatively." (See ref. 3, p. 249.) The methods used required concurrent compila- 
tion (multiprocessing) of both problems (ref. 3, p. 250). The next logical step in software 
development, which is a necessity for a general-purpose simulation laboratory, is the 
ability to process and/or compile more than one problem independently (multiprograming) 
on one general-purpose computer. 

With the advent of the next generation of digital computers, one envisions the utility 
of multiprograming in an effort to reduce cost in a general-purpose simulation laboratory. 
The removal of the restriction that the computer be dedicated to a fixed number of prob- 
lems imposes greater demands on the simulation programer. It is no longer sufficient 
merely to generate a working program, but it now becomes imperative that the most 
efficient way to mechanize a given simulation problem be determined. 

In transforming a simulation problem to a digital program, attention must be given 
to algorithms to approximate required mathematical functions. The analysis and devel- 
opment of integration algorithms or integration schemes shall be the body of this report. 

In recent years, there has been a general transfer of large simulation problems either 
whole or in part to digital computing equipment. As is true in any new field, the points of 
view usually cover a wide spectrum and the use of integration schemes in real-time sim- 
ulation is no exception as witnessed in references 1, 2, and 4 to 9. 

The general problem of differential-equation solving has led to many varied types 
of integration techniques which are usually complicated and often time consuming, the 
attempt usually being made to solve the widest possible class of problems most of the 
time. However, these techniques are often too slow and complicated for effective use in 
real-time problems. The approach taken in this paper is to consider only a rather 
restricted class of differential equations from a real-time point of view. The hope is to 
produce a fresh look at real-time simulation that is on a sound analytical and experimen- 
tal basis. 

Consideration is given in this report to the simulation of nonaerodynamic rigid 
vehicles with on-off control. The results however will in part be applicable to most 
flight-vehicle simulations. The accuracy requirement is to meet or exceed that of a 
corresponding analog solution (usually quoted on the order of 1 percent but not always 
acquired). The purpose of this report is to establish guidelines for meeting accuracy 
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requirements of this class of real-time simulation problems and to maintain a high 
degree of computer efficiency. 

The basic equations used to simulate the Gemini-Agena configuration are given in 
an appendix by Roland L. Bowles of the Langley Research Center. 

SYMBOLS 

The International System of Units (SI) are employed in this report. For conversion 
to U.S. Customary Units, consult reference 10. All systems used throughout this paper 
are right-hand, orthogonal, Cartesian coordinate systems. (See fig. 1.) 

A transformation matrix relating observer (Gemini) body axes and target 

(Agena) local vertical axes, dimensionless 

A column vector with components a, b, c, and d, dimensionless 

a,b,c,d quaternion elements (Euler parameters), dimensionless 

a^j matrix elements of A (i,j = 1,2,3), dimensionless 

B transformation matrix relating observer vehicle body axes to inertial frame, 

dimensionless 

B f transformation matrix relating target vehicle body axes to inertial frame, 

dimensionless 

b^j matrix elements of B, dimensionless 

bjj T matrix elements of B T , dimensionless 

D matrix relating the target vehicle axes to a line-of- sight axis system, 

dimensionless 

d ij matrix elements of D, dimensionless 

E enlargement or displacement operator, dimensionless 

e,f,g,h target quaternion elements, dimensionless 
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F matrix relating observer vehicle axes to target vehicle axes, dimensionless 

F 0 total initial fuel, kilograms 

F trans’ F att percent fuel used by translational and attitude control systems, 

respectively, percent 

fjj elements of the matrix F, dimensionless 

G(t) Green's function of linear gyroscope, dimensionless 

G^(t) associated Green's function of linear gyroscope, dimensionless 

g e acceleration at surface of earth due to gravitational attraction, 

9.84 meter s/second^ 

g D gravity acceleration at r s = r 0 , meter s/second^ 

H target vehicle rotational angular momentum vector, kilograms-meterVsecond 

Hi,H 2 ,H 3 components of target angular momentum vector H resolved into 

target body axes, kilogram -meter2/second 

h integration interval size, seconds 

I inertial matrix of observer vehicle, kilogram-meter2 

I"1 inverse of observer vehicle's inertia matrix, (kilogram -meter2)"^ 

Ijj elements of matrix I, kilogram -meter 2 

Ijj' elements of matrix I“l, (kilogram -meter2)~l 

I g p specific impulse of observer-vehicle reaction control system, second 

i imaginary unit, fi 

J inertia matrix of target vehicle, kilogram-meter^ 
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J“1 inverse of target vehicle’s inertia matrix, (kilogram-meter2)“l 

J ij elements of the matrix J, kilogram-meter2 

Jjj' elements of the matrix J"l, (kilogram-meter2)“l 

Ki thrust misalinements (i = 1,2,. . .,5), dimensionless 

Kq + 1 normalized target orbital angular momentum per unit mass, dimensionless 

gain constant used in maintaining orthogonality of transformation matrix 

L observer vehicle rotational angular momentum vector, 

kilogram-meter2/second 

Li,L 2 ,L 3 components of observer-vehicle rotational angular momentum L 

resolved into observer body axes, kilogram-meter2/second 

l column vector of partitioned direction cosine matrix, dimensionless 

Z,m,n partitioned components of direction cosine rate equations (components of l), 

dimensionless 


M moment vector acting on observer vehicle due to translational and attitude 

control reaction jets, joules 


Mi,M 2 , m 3 


components of vector M, joules 


X-components of attitude control moments, joules 


X- components of moment resulting from transla- 
tional control forces, joules 


( m 2)t 0± y “ co 

( M 2) T x± > ( m 2) T y± >( M 2 )t 


Y-components of attitude control moments, joules 


Y-components of moment resulting from transla- 
tional control forces, joules 
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II 


( M 3) T x± ’( M 3)t y± ’( M 3) T 


Z -components of attitude control moments, joules 


Z- components of moment resulting from transla- 
tional control forces, joules 


m 


mass of observer vehicle, kilograms 


^ m att>^ m tran 


mass changes resulting from translational and attitude reaction 
control jets, respectively, kilograms 


N 


total rotational angular momentum of observer vehicle squared, 
(kilogram-meter2/second)2 


N b 


number of significant binary bits, dimensionless 


n,N 


integers, dimensionless 


0(h n ) 


remainder in a series with lowest order term of h n 


P norm squared of A, dimensionless 

p,q,r components of observer angular velocity vector co resolved into observer 

vehicle axes, radians/second 

p t ,q|.,r^ components of target angular velocity vector a5^ resolved into target vehi- 
cle axes, radians/second 

p ff ,q ,r angular velocity components p, q, r measured by gyroscopes, 
radians/ second 


Q 

R 


norm squared of f, dimensionless 

relative range vector directed from center of gravity of target vehicle to 
center of gravity of observer vehicle, meters 


R 

Rios 


relative range rate, meters/second 


line-of-sight range vector, meters 
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Rg 0 components of relative range vector R resolved in observer 

vehicle axes, meters 

mean radius of earth, meters 

radius of target reference orbit, meters 

position vector of target vehicle in a coordinate system with origin at earth's 
center and coordinate directions always parallel to an inertial frame, 
meters 

polar coordinates locating the center of gravity of target vehicle with respect 
to inertial space (for elliptic orbits 6s = tt/2 at perigee), meters and 
radians, respectively 

skew-symmetric matrix that relates time rate of change of components of a 
constant magnitude quaternion to its components referenced to a rotating 
frame, radians/second 

Laplace transform variable 


components of observer vehicle thrust forces resolved in target local 
vertical axes, newtons 



magnitudes of on-off reaction control jets, newtons 


jTg k total forces along observer vehicle axes resulting from transla- 

tion and attitude control jets, newtons 


time, seconds 

target tangential velocity component (V = r s 0 s ), meters/second 

characteristic velocity of a vehicle in a circular orbit at radial distance r Q , 
meters/second 

total velocity change using translational control, meter s/second 


Ill 

AV^AVyjAVg components of velocity change using translational control along 

observer's X, Y, and Z axes, respectively, meters/second 

X,Y,Z components of relative range vector R referenced to a rotating coordinate 
frame with origin located at center of gravity of target vehicle, meters 

Xp,Yp,Zp coordinates of pilot's eye with respect to observer vehicle axes, meters 

AX,AY,AZ components of R los resolved in observer vehicle axes, meters 

X+,X- translation forward, translation aft, respectively 

Y+,Y- translation right, translation left, respectively 

Z+,Z- translation down, translation up, respectively 

x,r dummy variables of integration 

z z -transform variable 

a, j 3 azimuth and elevation angles of target with respect to observer vehicle axes, 

radians 

A forward finite -difference operator, dimensionless 

V backward finite-difference operator, dimensionless 

6 velocity perturbation variable, dimensionless 

e global error in angular momentum squared using Euler integration, dimension- 

less; also used as orthogonality error, dimensionless 

e r local relative roundoff error 

e tr local relative truncation error 

gyroscope damping ratio, dimensionless 

77 angle between l and o 5 , radians 
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e+,e- 


A 

\,r 




P 


0 + , 0 - 


X 


*,e,<P 


'P+,'P- 


’/ / C>^Cj0C 

ft 


w 

U) c 

w g 

Wo 

CO 


upward and downward pitch, respectively 

variable associated with gyroscope equation, radians/second 

target azimuth and elevation angle with respect to local vertical axis system, 
radians 

Euler angles which specify orientation of scaled target model axes with 
respect to simulation laboratory analog of line-of-sight axes, radians 

displacement perturbation variable, dimensionless 

right and left roll, respectively 

angle between l o and L, radians 

Euler angles relating observer vehicle axes to target-centered local vertical 
system, radians 

right and left yaw, respectively 

attitude angle errors used in attitude control system, radians 

Euler angles relating observer vehicle axes to an inertial frame, radians 

skew-symmetric matrix which relates components of a vector to their 
derivatives due to referencing a rotating reference frame (equivalent to 
ceX operation), radians/second 

magnitude of a), radians/second 

general characteristic frequency, radians/second 

gyroscope natural frequency, radians/second 

orbital angular rate for target reference orbit, radians/second 

angular velocity vector of observer vehicle with respect to inertial frame, 
radians/second 
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(•) 

On 

() T 

[4 

[-4 

( )(j) 


angular velocity vector of target vehicle with respect to inertial frame, 
radians/second 

denotes time derivative, second" * 
evaluation at time t = nh, dimensionless 
transpose operation, dimensionless 

transformation matrix of a simple rotation through any angle £ about X-axis, 
dimensionless 

the inverse and/or transpose of dimensionless 

jth derivative, second”^; also jth difference 


Axis systems: 


Xi,Yi,Zi an "inertial" system (origin at center of the earth X^Zi-plane defines orbit 
plane of target vehicle, passing through perifocus and Yi in direc- 
tion of orbital angular velocity) 


X,Y,Z target local vertical system (a rotating system, origin at center of gravity 
of target, Z-axis through center of earth, XZ-plane lies in target orbit 
plane. X-axis points generally backward in target orbit) 


X T,b’ Y T,b> z T,b 


x o,b’ Y o,b> z o,b 


x los’ Y los> z los 


principal axes through center of gravity of target (x T ^ 
forward, Z T ^ down, Y-j- ^ completes the right-hand 
system, docking ring aft) 

observer body-axis system (origin at center of gravity of the 
observer vehicle; Xq^ out nose of vehicle (forward), Y 0) k 
out right wing, Z Qj b (down) completes right-hand system) 

line-of-sight axis system (origin coincident with origin of 
observer body axes and X^ os passes through origin of 
target local vertical, Y^ os lies in the X 0; kY 0jb -pla.ne) 
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PROBLEM DESCRIPTION 


Background Discussion 

The physical systems considered in the analysis are those of a class of problems 
rather than of any specific problem in order to approach the widest possible degree of 
applicability. The physical systems of concern are those of manned and unmanned flight 
vehicles. The greatest variability of vehicle simulations lies in three parts: structural 

effects, the coordinates used for kinematical variables, and the type of control forces. 

The vehicles are considered as structurally rigid except in an example problem of a 
system of two elastically coupled rigid vehicles studied in the latter part of this report. 
The rotational kinematic variables chosen are, in most cases, relatively uniform, a wide 
variety being used in the translational kinematics. The direction cosines and the Euler 
parameters (unit quaternion) are considered as representations of rotations since they 
overlap most present-day simulations. The angular momentum components are used for 
the moment description, and they reduce to the angular rate equations for constant iner- 
tias. The translational variables are considered in the analysis only in very general 
terms. The applicability (or lack thereof) of the analysis of the translational motion to a 
specific problem is clarified in the section "Analysis." All statements apply to any 
flight-rigid vehicle simulation as long as the control forces are zero. The only control 
forces considered are derived from a reaction — jet control system. The complete 
analysis is then only for the real-time simulation of rigid vehicles with on-off control. 
However, inferences to other problems are made where possible. 

In addition to the on-off control forces, the effects of a gravitational field are also 
considered. The class of orbital-rigid vehicles with on-off control can be characterized 
by the frequencies associated with the motion. These characteristics are low frequencies 
(~10“3 radian/sec) associated with the steady-state trajectory motion, moderate frequen- 
cies («1 radian/sec) in attitude maneuvers, and usually high frequencies (~10 to 
1C)2 radians/sec) in the attitude-control system. After the analysis, two examples are 
taken from this class. The first is a Gemini-Agena rendezvous and station-keeping prob- 
lem (ref. 11), the mathematical model of which is in appendix A. The second is a varia- 
tion of this problem which includes an elastic tether (ref. 12). Also discussed are recent 
applications of the results of this analysis to several large simulations, some of which 
include aerodynamics. 


Fundamental Differential Equations 

Associated with every manned vehicle is a set of fundamental differential equations, 
the solutions of which yield the vehicle’s orientation, angular rates, position, and veloci- 
ties. The description of differing vehicles mainly changes the coefficients and the inho- 
mogeneous term (including the control system) of these fundamental equations. 
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Axis transformations. - Two methods of computing the transformation matrix are 
commonly used, both of which are considered in the analysis. These methods use the 
direction cosine rate equations and the Euler parameter rate equations. 


The first and probably the most common method is use of the direction cosine rate 
equations. The nine direction cosine rate equations are partitioned into three related 
sets all of the form 


— — 


— ~ 


— 

i 


0 r -q 


l 

m 

= 

-r 0 p 


m 

n 


f 

1 

O 

l 


n 


(1) 


Each related set has the property of normality 

Q = + m2 + n2 = 1 


( 2 ) 


which is a constant of the motion. 

Equation (1) can be written as 


I=m 


(3) 


where l is the vector of components l, m, and n and is the skew-symmetric 
matrix in equation (1). 


The Euler parameters are obtained from the following matrix differential equation 
(eqs. (A35) to (A39) with Kg = 0) 


a 


0 r -q -p 


a 

b 

_ 1 

~ | 

1 

0 

1 

►a 


b 

c 

d 

2 

q -p 0 -r 

p q r 0 

_ 


c 

d 


Equation (4) can be written as 

1 = SA 


(4) 


(5) 


where A is the vector (a, b, c, and d) and S is the skew-symmetric matrix in 
equation (4) including the factor The orthonormality condition is (eq. (A3 9) for e = 0) 
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P = a 2 + b 2 + c 2 + d 2 = 1 (6) 

Hence, P given by equation (6) is a constant of motion. 

Attitude co ntrol system .- The simulation of the attitude control system is comprised 
of computing angles from the gyroscope outputs: 


<t>c = P g 
$c = *lg 

*C = r £ 


(7) 


where p g , qg, and rg are the measured body angular rates and are usually given by 
the approximate linear gyroscope equations: 


Pg + 2 ?g“gPg + w g 2 Pg = w g 2 P 
q g + 2CgO) g qg + w g 2 q g = w g 2 q > 
r g + 2?gCi> g r g + w g 2 r g = o) g 2 r 


( 8 ) 


where p, q, and r are the true body angular rates as calculated from the moment 
equations. 

Mo ment equations . - The moment or angular momentum rate equations are written 
as the following matrix equation (eqs. (A43) to (A45)) 


Ll 


0 r -q 


Ll 


Mi 

L2 

= 

-r 0 p 


L2 


m 2 

l 3 _ 


►Q 

O 

\ 


L3 


m 3 ^ 


where from equations (A40) to (A42), 


P 


111 

112 

113 

-1 

Ll 

q 

= 

*21 

l 22 

123 


l 2 

r 


hi 

*32 

X 33 


l 3 


(9) 


( 10 ) 
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Equations (9) and (10) can be written as 

L = OL + M (11) 

and 

uJ = I -1 L (12) 

— T 

Forming the scalar product with L and equation (11) yields 

— T— T — — T — 

L L = L fiL + L M (13) 

When M is null, the right-hand side of equation (13) vanishes since L is orthog- 
onal to £2L. This relation implies that N (the total angular momentum squared) given 
by equation (14) is a constant of the motion for M null: 

N = L T L (14) 

Translationa l equations.- The translational equations are found in appendix A 
(eqs. (A4) to (A10)). They are characterized by very low frequencies 

* 10“ 3 radian/sec) for nonthrusted motion. For thrusted motion, the frequencies 
become moderately high (^o c ~ 1 radian/sec from rotational motion coupling^ but have 
little effect on the solution character. ( The effect is a small-amplitude oscillation with 
(jOq ~ 1 radian/sec modulated with a large-amplitude oscillation with 
c Oc ~ 10~3 radian/sec). The constants of motion for no thrust would be the total energy 
and orbital angular momentum. 


ANALYSIS 

The errors considered in the analysis are of two types: the error introduced by 
using a finite-word structure (round-off error), and the error in using approximate for- 
mulas for limiting processes such as integration or differentiation (truncation error). 

The following sections begin with a discussion of local round-off and truncation error. 
These local errors are related to the accumulated or propagated error by error domi- 
nance. The concept of error dominance and its region applicability are discussed in the 
next section. Error dominance is then the guide to the analysis of the equations of motion. 
The analysis of each related set of equations is followed by a discussion of the results of 
experiment. Before beginning this task, it first seems advisable to establish the ground 
rules under which the analysis is performed. 
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Characteristic of any nontrivial simulation problem is the evaluation of many com- 
plicated algebraic expressions for visual display drives, readout, and the derivatives of 
the systems variables. The performance of the integration is, with few exceptions, the 
minor part of the computation (only a few arithmetic operations per integration). Hence, 
integration schemes which require the least number of derivative evaluations per inte- 
gration interval are desirable (that is, multistep methods are preferred over the single- 
step methods). Equally important is the need to minimize the number of integration steps 
needed to obtain the problem solution (thus, the arithmetic unit is used as little as pos- 
sible per problem second). Clearly, the interval size is to be made as large as possible 
in order to reduce the amount of computation required and small enough to provide 
"proper" response to the pilot input. Reference 13 indicates that intervals on the order 
of 100 milliseconds is adequate for sampling pilot input. However, the increased accu- 
racy for smaller intervals is not clear because of roundoff error in that study. (The 
maximum word length was 7 bits.) For the remainder of this report, the maximum 
interval is assumed to be 50 milliseconds or 20 samples per second which should be 
more than adequate. (This sample rate is a factor of 2 more often than that suggested by 
ref. 13.) The minimum interval size is limited by the roundoff error and is discussed 
later in this report. 

Two methods are used to study the propagated truncation error. The z-transforms 
(ref. 14) are used to study the linear portions of the problem. The changes in the con- 
stants of motion are also used to indicate the amount of error introduced in the problem 
solution by various integration schemes. 

Roundoff and Truncation 

Analytica l development .- The errors of numerical computation are of two classes: 
errors due to representing numbers by a finite word size (that is, roundoff) and errors 
due to neglecting higher order terms in approximate techniques (that is, truncation). 
Although the effects due to truncation can, more or less, be determined, the propagated 
error associated with round-off can seldom be rigorously studied for practical problems. 
Rather, the propagated roundoff error is studied through model building and experimenta- 
tion. (See ref. 15, p. 41; ref. 16, p. 305.) 

The equations of motion are approximated (in most cases) by finite difference equa- 
tions of at least the same order of the equations of motion for consistency. (See ref. 15, 
p. 224.) Any finite-difference approximation of order greater than the order of the equa- 
tions of motion introduces extraneous roots. (The language of linear equations is used 
here; the meaning should be clear.) In addition to the extraneous roots are other roots 
which approach the true roots of the differential equation in the limit as the interval 
approaches zero (that is, h — 0) so that nh is finite and equal to t. If the extraneous 
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roots are unstable (ref. 15, pp. 185, 242), the approximate solution is sensitive to small 
errors (for example, starting error or roundoff error) introduced in the calculation. 
(See ref. 17, p. 11.) This sensitivity (called numerical instability) usually occurs when 
the sample rate and the systems frequencies are of the same order of magnitude. 


Implicit in the following sections is the assumption of numerical stability from 
which follows the concept of error dominance. Error dominance can be stated as follows: 


The propagated error will be dominated by 

1 


J roundoff 


i 


numerical stability, the local 


local 


J truncation l 
\ roundoff j 


J roundoff 
|truncationJ 


error if, in the region of 
error is in magnitude much larger than the 


|^truncationj 


error. This definition leaves an uncertain region between these 


regions where local roundoff and truncation error are within an order of magnitude of 
each other. Hence, these local errors compete in forming the propagated error. 


In subsequent sections, error dominance is used to determine the error that is 
most likely to affect the solution. The equations of motion are then analyzed on this basis 
and compared with experimental results. However, the magnitude of the local relative 
roundoff and truncation errors must be determined first. 


The local roundoff error is subdivided into two components, as in reference 15 
(p. 36), the inherent error and the induced error. Extensive discussions of these errors 
are given in reference 15 and are not repeated here. Instead, a few simple notions about 
roundoff error are given and the reference is cited for a more extensive discussion. For 
the purpose at hand, the inherent error can be considered as the error in evaluating the 
derivatives of the system. Upon integration these errors are reduced in magnitude by 
the integration interval h. The induced error is introduced through the actual perfor- 
mance of the integration and is not reduced by h but is always of the order of the least 
significant bit of the finite word structure (ref. 15, p. 35). 

Before deriving an approximate model for the local relative roundoff and truncation 
error, consider a few simple operations in a hypothetical normalized 9-digit floating- 
point decimal machine. The word structure for a number is represented schematically as 


A = (.aaaaaaa. . . .) x 10 


exp a 


0.1 = A X 10“ EXP A < 1 


where the a terms are digits of values 0 to 9 and EXP A is the integer exponent of A. 
The decimal machine word structure is represented as 
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exp a 


.aaaaaaaaa 


for a 9 significant digit machine. Note that the relative error of this machine represen- 
tation of A is bounded by 10“ ® (the relative value of the first significant digit being 
neglected by the decimal machine). The inherent error can largely be controlled by care- 
ful programing. Taking differences of nearly equal numbers results in a large inherent 
error, which is demonstrated schematically as follows: 


EXP 

.XXXXYYYYY 


EXP 

.XXXXZZZZZ 


A - B = EXP - 4 .DDDDDRRRR 


Roundoff 


D = Y - Z 

The roundoff error represented by the R digits move into the register when left shift 
is performed after subtraction. Thus, if this operation or similar ones are used, the 
inherent error may be significant. With careful programing, the inherent error can be 
neglected because of reduction by the interval h upon integration. 

If it is assumed that the inherent error is controlled, the induced error must be 
studied. The source of the induced error resulting from a simple numerical integration 
procedure is shown as follows: 
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Shift right on AX n exponent until equal to the X n exponent. Then add: 



No. of zeroes = Inherent error 

EXP X - EXP X - EXP h 

exp x I xxxxxxx N x N x N 


X N = X + X A 

The variables X^, Xi, and h are shown (with the inherent error in &) in their decimal 

word representation. (h = 10 EXP h is assumed for simplicity.) The integration is per- 
formed by (Euler integration is used for simplicity) 

X n+ i = X n + hX n 

X n+ 1 = X n + AX n 

as shown schematically. The AX is computed, and is followed by a shift right (usually) 
until the exponents of X and AX are the same. The addition is then performed. Note 
that the induced error is that part of the AX word shifted out of the register. Thus, 
the relative weight of the induced error in X is always 10“^ (the relative weight of the 
rounded digit). 

The control of the induced error could be performed by using a longer register and 
partial double precision as follows: 



h = EXP h + 1 . 100000000 
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AX n = hX n = EXPx + EXP h .X A X A X A X A X A X A X A R A R A 


Shift right on AX n exponent until equal to the X n exponent and then add 


X n + AX n = 


EXP 


X 


.xxxxxxxxx 


xxxxxxxxx 


Induced error 


EXP X 


.000000X A X A X A 


X A X A X A X A R A R A 


RRR 


Inherent 

error 


EXP-< 


. xxxxxxx N x N x N 




X N = X + X A 


The inherent error is 


Rjvj = X + R a 


The induced error is 


Rj^ — X + R 

The register containing X, if sufficiently increased in length, would contain mainly the 
inherent error reduced by the integration interval size h. The derivatives X may 
still be evaluated by using the normal word size and two words for X; this mode of 
operation (partial double precision, ref. 15, p. 94) normally reduces roundoff error sig- 
nificantly. Examples are shown later in the text. 

In order to make effective use of error dominance, the relative magnitudes of the 
local roundoff and truncation errors are needed. It will henceforth be assumed that good 
programing practices are followed and the inherent error can be neglected. A simple 
model is derived for the local relative roundoff and truncation error (inherent error being 
neglected) for second-order integration techniques and a normalized N-bit floating-point 
machine (binary). 


19 





It is assumed that the jth derivative of the state variables of the system are given 
by the relation 


X^ = w^X 


(15) 


where the characteristic frequency a> c is, in general, complex. 

Consider a Taylor series expansion of X from the point t = rih to t = (n + l)h 


C n + 1 = X n + h *n + \ h 2 X n + 1 h 3 X n + . 


Substituting equation (15) into equation (16) results in 

Xn +1 = (l + ho> c + 1 h 2 0J c 2 + | h 3 w c 3 + 


n+1 - (l + + x h^ w c^)^r 


The solution using second-order integration is 

X 

and the local truncation error T is 

-itAa S 


T * | h'V'Xn + o(h 4 ) 


The change in X to second order relative to X is 


AX n = hco c M l + ho)c| 


X, 


■n 


(16) 


(17) 


(18) 


(19) 


( 20 ) 


which is approximately the change at the nth step of the normalized portion of the com- 

rvp 

puter word representing X since X n ~ 2 x . Because of the finite word structure 
of the machine, the relative change in X to second order is approximated in the calcula- 
tions by 


f AX- 


X*T/a * hC ° c V + ^ ha>< y ±2 Nb 1 


(21) 


where is the number of significant binary bits. Hence, the local roundoff error R 

is bounded by 


-2 -Nb " 1 |x n | ^ R 22" Nb " 1 |x n 


( 22 ) 
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The magnitude of the local relative roundoff error e r is defined as 



which has a maximum of 


( 23 ) 


,-Nb-l 


e r = 


ho> r 1 1 + — hu)f 


(24) 


and is a measure of the greatest significant change in X 
deleted from the calculation due to rounding. 

The magnitude of the local relative truncation error 
being neglected, is 


that can occur and still be 


e tr , the remainder 



e tr = 


_T _ 
AX^ 



(25) 


The two local relative errors given by equations (24) and (25) are the limits of 
resolution imposed on X by the word structure and integration scheme, respectively, 
as seen by the variable X. 

In closing this section, it seems instructive to indicate the effects of statistical 
modeling on the local relative roundoff error. Two commonly used models for the dis- 
tribution of roundoff error is to assume that roundoff is uniformly distributed or Gaussian 
distributed. Implicitly assumed is that roundoff error at the nth step is independent of 
the error at the previous steps. The uniform distribution is over the interval (0,er) and 
the 50 percentile occurs at the midpoint, namely l/2e r . That is, 50 percent of the errors 
will be less than l/2e r for the assumed uniform distribution. The 50 percentile for the 
Gaussian distribution occurs at 0.68STD (STD is the standard deviation) where e r = 4STD 
is assumed. The 50 percentile point is then 0.17e r . 

Discussion of app lications .- The models for the local relative roundoff and trunca- 
tion error were derived by assuming the state variables of the system and their deriva- 
tives were linearly related. This assumption implies that the system has a degree of 
smoothness in its motion so that over short intervals of time (on the order of h), the 
motion is quadratic. The degree of variation of the characteristic frequency as a function 
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of the state variables and time clearly determines the validity of this local model. Dis- 
continuities that arise from thrust destroys this linear relationship and this difficulty 
will be dealt with later in the section "Moment Equations." 

The approximate local relative roundoff (the upper bound e r ; the uniform 50 per- 
centile e r U; and the Gaussian 50 percentile e r G) and truncation errors are shown in 
figure 2 for = 27 bits. These local errors are estimates of magnitudes only and no 
information is supplied for additivity (which is machine and program dependent). 

The manner in which these errors propagate is a more complex problem than this 
simple analysis can ascertain. However, this analysis does indicate that the error most 
likely to enter the problem (error dominance) and an estimate of an upper bound of the 
propagated error can be found. (See ref. 15, p. 36.) 

As seen in figure 2 with = 27 bits, the portion of the computation including the 
attitude -control system (hu> c ~ 1^ and moment equations and axis transformations 

(hu> c ~ 10" 2 ) could be greatly affected by truncation with little concern for roundoff. If 

low-order integration schemes are to be used in these portions of the simulation, great 
care must be employed to determine any possible adverse effects. However, the trans- 
lational equations (ho> c ~ 10 - ^ are relatively insensitive to truncation error but roundoff 
error could be a prime difficulty. 


Axis Transformations 

Analytical dev elop ment .- The transformation matrix is computed either directly by 
integrating the direction cosine rate equations or by the Euler parameter rate equations, 
both of which must be subject to the restrictions of orthonormality. By using the steepest 
descents, the rate equations could be augmented with constraints to maintain orthonor- 
mality. (For example, see eqs. (A35) to (A39).) Although orthonormality would be main- 
tained, phasing errors could become appreciable. In this vein, an integration scheme is 
chosen that closely approximates the orthonormality condition. This scheme will then be 
used in conjunction with a normalization procedure to maintain orthogonality of the trans- 
formation matrix. 

First consider Euler integration to solve the direction cosine rate equations given 
by equation (3) 


^n+1 _ ^n + hZ n - ^n + h^n^n 


(26) 


The norm squared (eq. (2)) by using equation (26) becomes (where = coXZ) 


22 



(27) 


Qn+1 _ Qn + (^ n ^n) (^n^n) - Qn + (hw n sin rj n ) Q n 

where rj n is the angle between the l n and the 5J n at t = nh. From equation (27), 
the per step error due to Euler integration is 


AQ n = (hw n sin 77 n ) Q n 

Similarly, the Adams-Bashforth integration formula yields 

^n+1 = ^n + 2 ^n^n " ^n-l^n-l) 


(28) 


(29) 


and the per step error of the norm is 


where 


AQn — hf n fin-l^n-1 + U U 


U - 3f2 n Z n - ^n-l^n-1 


(30) 


A striking difference in the change of the norm is apparent for the two schemes, 
as seen in equations (28) and (30). Euler integration gives rise to a positive definite 
change in the norm as opposed to Adams-Bashforth integration which does not have a 
definite norm change. Hence, Euler integration always increases the norm for each 
integration step. For example, consider the ratio of the change in the norm squared to 
the norm squared with Euler integration for hu> ~ 10"^ radian 

■^^2 «h^u>2 » 10“^ = 0.01 percent 

Qn 

Hence, an error of 0.01 percent is committed per integration step. Clearly, a solution 
for a reasonable length of time is not acceptable even for moderate frequencies 
(u *1 radian/sec) and interval sizes (h * 10"^ sec). 

Consider equation (27) when 53 is a constant. The z-transform of equation (27) is 
(sin t) is also constant to a very good approximation) 

(z - 1 - h^w^sin^T^QU) = 0 (31) 
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II II llll 


II II 


l ii ii mi 


which has one root 

z = 1 + h 2 o> 2 sin 2 77 (32) 

The solution to equation (27) is then 

Q n = z n = e 11 log e< z ) (33) 

By using the assumptions that 

h 2 a> 2 sin 2 T} « 1 (34) 


and 


log e (l + X) ~X-|x 2 +|x 3 -ix 4 (35) 

an approximate solution (the first nonvanishing term in approximation (35) being retained) 
for the norm squared is 


Q * e 


h(cu sin 77 ) z t 


(36) 


when Euler integration is used. 

The solution of the direction cosine rate equations for Euler integration is found by 
taking the z-transform of equation (26) with a constant. The z-transformed equation 
can be written as 


[(z - 1)1 - hfi]z(z) = 0 

where the identity matrix I in this and similar equations is suppressed in the sequel. 
Hence, the equation is written 

[(z - 1) - hn]f( Z ) = 0 (37) 

The characteristic equation of equation (37) is 

(z - l)[(z - l ) 2 + h 2 u> 2 ] = 0 (38) 
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for which 


Z 1 = 
z 2,3 



are solutions. 

The solution to equation (26) is then 

l n = ci + c 2 e n lo ged+^h) + g^n log e (l-icoh) 


( 39 ) 


(40) 


where the Cj terms are constant vectors that depend on initial conditions. Again with 

approximation (35) and by assuming that | ihaj J « 1, equation (40) is rewritten by 
retaining the first nonvanishing terms in damping and frequency errors as 




l- Ci + e 


c 2 e 


„ h 2 a,3) t 

CO / t 


tl 


+ c 3 e 


(41) 


Equation (41) is in agreement with equation (36). 

The solution for the direction cosines using Adams-Bashforth second-order inte- 
gration is found by taking the z -transform of equation (29) for constant £2 (where the 
identity matrix is implied): 


z(z - 1) - § £2(3z - 1) 


The characteristic equation of equation (42) is 


i(z) = o 


z(z - l)||z(z - !)] 2 w 2 (3z - 1)^ 


= 0 


Two roots of equation (43) are by inspection: 


zi = 0 


Z 2 = 1 


and the four other roots are given by 


(42) 


(43) 


(44) 

(45) 


|z(z - l)j^ = -c^(3z - 1)' 


(46) 


25 



where 



Equation (46) can be reduced by taking the square root of each side: 


z(z - 1) = ±ic(3z - 1) (47) 

The roots of equation (47) are then given by 



The four roots represented by equation (48) are two sets of complex conjugate pairs. 
This fact can be verified by substituting -i for i. 

The real and imaginary parts of the roots given by equation (48) are found by 
defining 


ue* a = (l - 9c 2 ) + i2c 


It is seen from this definition that 


where 


u l/2 e i (: ’'/2 - J^i _ 9c 2 ) + i2cj 

u^/^e”* 0 /^ = j^i _ g c 2) _ j2cj 

u = (l - 14c 2 + 81c 4 ) 1 / 2 


a = tan 


Hence, the roots given by equation (48) are 


-l/_2c 


VI - 9c' 


5 3,4 = 2 


Z 5,6 = o 


^1 + u V 2 cos ± i ^3c + u 4 / 2 sin ^ 

1 - u 4 / 2 cos ± i ^3c - u 4 / 2 sin ^ 


(49) 

(50) 


( 51 ) 


J 
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If c « 1 is assumed ^for c = » 10’^ radianj , equations (49) and (50) can, terms 

in c5 and higher being neglected, be approximated by 

u = 1 - 7c2 + 16c^ (52) 


a = 2c ^1 + -^ C 2j 

and similarly 

u 1/2 = i_ 7 c2 + 15 c 4 
2 8 

Sin f =C ( 1+ f c2 ) 


(53) 


(54) 

(55) 


£= i _I C 2 - 61 c 4 

2 2 


(56) 


Substituting equations (54), (55), and (56) into equations (51) results in 

z 3 4 * 1 “ 2c2 " 2c^ ± i(2c + 2c ^ 
z 5 6 * + ± i( c “ 2c^ 


(57) 

(58) 


The two roots given by equation (58) are quickly damped for ha> small or equivalently 
for c small. If the heavily damped roots are neglected and approximation (35) is used, 
the solution of the direction cosines for the Adams-Bashforth second-order integration 
scheme is given as 


+ e 


h 3 o> 4 , 

16 


i (oo+— h2u>3) t -i (w+— h 3 co3 

l 12 / _ V 12 / 

c 2 e c 3 e 


(59) 


where Cj, C 2 , and c 3 are the same as in equation (41). By comparing the solution 
given by equation (59) to that given by equation (41), the solution for the norm squared 
for Adams-Bashforth integration is by analogy with equation (36) 


Q 


8 

e 


h3q)4 

' 16 


(sin rj)^t 


( 60 ) 


when second-order Adams-Bashforth integration is used. 
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In order to make comparisons of the direction cosine formulation with the Euler 
parameter formulation, the Euler parameter rate equations will now be solved under the 
two integration schemes. Note that the z-transformed characteristic equations (eqs. (32) 
and (43)) are often similar in form so that the factors already obtained can be used in the 
following analysis and will reappear in later sections. 

If Euler integration is applied to equations (5) and (6), 


■A-n+1 ~ An + hS n A n 


2 ,., 2 


AP n = 


h o> 


(61) 

(62) 


From equation (62) for co constant, 


h^t 


P ~ e 


which is the same result as that obtained from the solution of equation (61) 

. ( x) h^o;2\ 


ho>* 

A = e 8 






2 24 / ^2 24 

a^e + s^e 


(63) 


(64) 


where again approximation (35) has been used and only the lowest nonvanishing order in 
amplitude and frequency errors has been retained. The a^ terms are constant vectors 
that depend on the initial conditions. 

Similarly, for Adams-Bashforth second-order integration, 

A n+ i = + T^(3SnA.n - S n _]A n _i) (65) 


2 2 

AP n = "hAn^Sn.iAn.i +--^3S n A n - S n _]A n _i^ 


( 66 ) 


Again, the solution of equation (65) is found by z-transforms for co a constant (that is, 
S is constant). 


z(z - 1) - | S(3z - 


1 


A(z) - 0 


(67) 
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1 _ . . 

where the identity matrix is implied. By letting c = — , the characteristic equation of 
equation (67) can be written as 


[z%(z - l)^ + c2(3z 



( 68 ) 


which is of the same form as equation (46). Hence, the roots are of the form of equa- 
tions (57) and (58) where again the roots of equation (58) are heavily damped. The solu- 
tion of equation (65), the heavily damped roots being neglected, is approximated by 


A = e 


4 h^w^tj 

256 


i^+-^h2ct>3]t 

2 96 


aie 


+ a 2 e 


-i(w+JLh 2 o> 3 

2 96 


(69) 


where the a.j terms are the same as those in equation (64). From equation (69), the 
solution for the norm squared for Adams-Bashforth second-order integration (eq. (66)) 
is found to be 

8 

P «e 



Discu ssion of experiment .- The results of this analysis are summarized in table I. 
The growth time constants and frequencies for both the direction cosines and Euler 
parameters when Euler and Adams-Bashforth second-order integration are used are com- 
pared with the parameters of the corresponding differential equations. 

The effects of the truncation error (hence, the propagated truncation error) for the 
direction cosine and Euler parameter rate equations has been determined for constant aJ. 
In all cases it is seen that frequency of the system is only slightly affected. The most 
prominent behavior (for oJ constant) is seen to be the amplitude growth which is 
reflected in the norm. 

The solutions to the norm of the direction cosine and Euler parameter rate equa- 
tions using Euler and Adams-Bashforth integration were obtained both analytically and 
on the computer. Euler integration was used to start the Adams-Bashforth scheme. The 
solutions were obtained for 


P = q = 


1 radians 
\J 2 second 


sin 77 = 1 
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initially and 


h = second 

16 

These solutions were compared with the true norm which is unity and the resultant errors 
are shown in figures 3 to 6. As seen in these figures, the solution parameters as given 
by table I are accurate. (Also shown are the starting errors caused by the Euler inte- 
gration starting formula.) 

The solution parameters in table I clearly show that the Adams-Bashforth second- 
order integration and the Euler parameters are superior to the other possible combina- 
tions. The time constant is better by a factor of 16 and the frequency error is also the 
smallest. This combination is sufficient for many simulation problems as witnessed also 
in figure 6. ^With h = 2~4 sec, j’cU | = 1 radian/sec for intervals of 100 seconds is 
more accurate than is required for most problems. j However, note that starting errors 
can be appreciable for nonzero p, q, and r when Euler integration is used to compute 
starting values. 

The assumption that u> is constant used in the analysis is not realistic. In gen- 
eral, aJ is a function of time with magnitude varying between zero and some co max . 

When w is small, roundoff error will again affect the solution. The magnitude | cij| 
is small during coasting periods of the flight where the control system is normally in 
attitude hold and the vehicle exhibits limit cycle motion inside a deadband. During this 
type of motion, roundoff error is the dominant error entering the solution. The main 
concern is placed on the stability of the limit cycle and usually the amount of fuel used 
during this portion of the flight. For most manned vehicles, the limit cycle motion is 
stable for a = 27 bit machine for h = 15 msec (this value is not a known absolute 
lower bound) and the experimental results are presented later. 

The w is, in general, time varying and the analysis is not applicable. Euler 
integration has the property that the change in the norm is always positive definite (see 
eq. (28)) and is always unstable. The Adams-Bashforth second-order integration is not 
positive definite and the norm may actually decrease for time varying u\ (This effect 
is seen later in another section.) 

Typically in flight simulation, attitude maneuvers are performed by sustaining some 
constant angular rate over a maximum time interval. An estimate of the adequacy of the 
two formulations can then be made by using the parameters in table I. One then has a 
basis for choosing a formulation and integration scheme for a particular application. 


30 



Attitude Control System 


Analytical development .- The response of the attitude-hold mode is strongly 
dependent on the quality of the solutions of equations (7) and (8), the attitude angle error, 
and the gyroscope outputs. If it is assumed that the solutions of equations (8) (the gyro- 
scope outputs) are truly representative of the angular rates of the vehicles, conclusions 
can be made as to the integration scheme to be used for equations (7). Consideration is 
only given to the first of equations (7) 


cp c = p. 


g 


(71) 


since the analysis holds for the 6 C and i// c equations as well. 

For an on-off attitude control moment, the angular rate is approximated by 


J g 


“TTT ^ 1 


dt 


(72) 


The moment of inertia is a slowly varying function of time; hence, over short-time 
intervals it may be considered to be constant. Since Mj is either zero or a constant, 
the rate given by equation (72) is a constant or linear function of time, respectively. 
Visually, 


/Mi 

p g = P g ,o + [I* 


(73) 


where the t coefficient is zero or a constant value and the subscript o denotes initial 
conditions. The solution to equation (71) is either a linear = 0^ or quadratic 

^ 0) function of time, as 


1 M ll *2 


*c = <t>c,o+P g , 0 t + 2i-)t 


(74) 


Consider Adams-Bashforth Lth-order integration as applied to equation (71) 


* c n+l 


c + h. 


L-l 


X r i V]p g,n 

1=0 


(75) 


where the yj values are given in reference 15 (pp. 192-193). Since Pg is at most a 
linear function of t, L = 2 in equation (75) is an exact representation of <p c under 
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these assumptions (eqs. (73) and (74)). This relationship is true since all higher orders 
of Pg vanish (ref. 18, p. 9). 

The equation for Pg is considered to be representative of equations (8). 

Pg + 2 ?g w gPg + w g 2 Pg = w g 2 P (76) 

The forcing function in equation (76) is a relatively smooth function of time. Even so, 
with typical parameters ( 'u>g ~ 10 2 radians/sec and £g ^0.8^, the numerical solution of 
equation (76) is not well behaved for simple numerical techniques when a reasonable 
interval size (h ~ 10“ 2 sec^ is used. 

The solution to equation (76) has been the object of research for many investigators. 
Examples of these efforts are given in references 4 to 7. The results of several such 
studies are summarized in reference 5. In view of the conclusions of reference 5, the 
next matter to be considered in this paper is to derive a method by which equation (76) 
can be solved. The complete derivation is contained in appendix B and only the high- 
lights appear in this section. 

The solution to equation (76), as found in appendix B, is 


where 


p g (t) = j* G(t - t) p(t) dr 


w t 


G(t) =— JL e"S T sin Tr 


(77) 


(78) 


with 




1/2 

COg 


« = «g«g 


and G(t) is a solution of the equation 

G(t - r) + 2?gW g G(t - r) + cc g 2 G(t - r) = Wg 2 6(t - r) 

and 6(t - r) is the Dirac delta function. (See ref. 19, p. 255.) As shown in appendix B, 
the solution given by equation (77) is valid even when p is a nonlinear function of p . 

o 
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(p is a nonlinear function of Pg for most control systems.) For further discussions of 
nonlinearities, see reference 20. 

A step-by-step evaluation procedure for solving equation (77) is found to be 
(eqs. (B36)) 




■ 


— ' 

1 

P g (t + h) 

X 

1 

CD 

II 

cos Fh 

sin Fh 

Pg(t) + ^ 

^h 

G(-x)p(t + x)dx 

'o 

A(t + h) 


-sin Fh 

cos Fh 

A(t) + ^ 

h j. 

G'(-x)p(t + x)dx 

0 

_ 






. _ 



If p(t) is known on the interval (t, t + h), equations (79) could be evaluated exactly. In 
practice, p(t) is known to the accuracy of the integration scheme used to obtain it. 
Assume that p(t) is known to within a linear function (Euler integration) on (t, t + h) 
(see eq. (73)) 

p(t + x) = a + bx (80) 

where 0 = x = h and a and b are constants over the interval that depend on p(t) 
and its derivatives. By substitution of equation (80) into equations (79), equations (79) 
can be written as 


p g (t + h) 

= e-^h 

cos rh 

sin Fh 

p g (t) + all + bl2 

A(t + h) 


-sin rh 

cos Fh 

A(t) + al^ + bl 2 t 


where 



(81) 


(82) 
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II 


Ill II I 


Evaluation of integrals (82) results in 

r 


~\ 


ii = -\ ?G(-h) + rG f (-h) 
w/ _ 


g 


- 1 


Il t =- L rK Gt( " h ) * rG( " h) 

w g L 


I 

r 


2 = — ^JhijG(-h) +hrG t (-h) - Hi - iV 

COg L 

2 f = ^-(h^G t (-h) - hrG(-h) - Hl f + 
w g L 




(83) 


The integrals (82) are constants that need only be evaluated once and then used as coeffi- 
cients in equations (81). For each interval step, only the coefficients a and b of 
equations (81) need to be evaluated. An obvious choice for a and b is 


a = p(t) 
b = P(t) 


(84) 


However, more accurate approximations for b can be made; for example, 

b = |[3p(t) - p(t - h^j (85) 

It is also a simple matter to make more accurate approximations to p(t) than that given 
by equation (80). 

When gyroscopes are involved, a more accurate approximation for <p c can be 
made than that given by equation (75). This approximation is made by integrating p g as 
given by equations (81) and expanding the results in a manner similar to this treatment. 

Discussion of experiment .- The linear gyroscope equations were solved by using 
the Green's function technique for a step and ramp input. The parameters were 
Wg = 125 radians/sec, = 0.8, and h = 1/64 sec. The response, as seen in figures 7 
and 8, were unchanged when the interval size was doubled. This result is expected since 
it is an analytical solution. Note that the solution even for discontinuities occurring at 
the beginning of the integration interval is still exact. 
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An independent study (ref. 8 ) has shown the h ~ 10“^ sec (10 000 solutions per 
second) is required for the Adams-Moulton second-order method to obtain a reasonable 
solution. This requirement is obviously a difficult one for simulation in real time. 

The Green's function technique employed to solve the gyroscope equations is exact 
when the forcing function is linear over the integration intervals. Higher order terms 
can be carried in approximating the forcing function with little added computer time. 

This method should be extended to include such effects as nonwindup limiting (piecewise 
linear gyroscope equations). This technique can also be extended to other linear and 
piecewise linear systems. 

Recently, Giese (ref. 21) has used a similar method for solving the linear part of 
control systems. He has developed an algorithm for computing an approximate Green's 
function for the general case. Although the approximate Green's function is step size 
dependent as can be seen in figure 12 of reference 21 , the algorithm approach is in the 
right direction. The Giese algorithm, however, does appear to have convergence problems 
as seen by comparing the results from second-order Adams-Bashforth solutions and the 
results of the Giese algorithm shown in figures 11 and 12, respectively, in reference 21. 


Moment Equations 

Analytical development .- The analysis of the moment equations is accomplished in 
two steps. First considered is the case of no external torques for which the total angular 
momentum squared (eq. (14)) is a constant of the motion. This analysis is followed by a 
discussion of applied on-off control torques. Throughout the analysis, it is assumed that 
the angular momentum and angular velocity are measured in the principal axes. The 
inertia matrix and its inverse are then diagonal. This procedure greatly simplifies the 
analysis without loss of generality. 


The moment equations are not linear. This nonlinearity is seen by writing the 
moment equations in component form, the angular velocity components being replaced by 
their functional form given by equation ( 12 ), 


-\ 


Li = rL 2 - qL >3 + Mj = L 2 L 3h^ - 


L 2 = pL 3 - rLi + M2 = L 3 L 1 M- - + M 2 } 


L 3 - qLj - pL 2 + M 3 = LiL 2 A-i -i-j + M 3 

1*22 hi 


( 86 ) 
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Although a solution to the general problem has not been found, two solutions have been 
approximated for two important types of vehicles. The first is vehicles with "near" 
cylindrical symmetry (that is, only two principal inertias are nearly equal), and the sec- 
ond is vehicles of "near" spherical symmetry (that is, all three inertias are nearly equal). 
For types of near cylindrical symmetry, the moment equations are approximately linear 
and z-transforms can be used. The moment equations are nonlinear for "near” spheri- 
cal symmetry but an approximate solution can be found for the Euler integration scheme. 


No external torques: The moment equations have the same form as the direction 
cosine rate equations for no external torques (that is, Mj = M2 = M3 = 0). It follows 
that the per step error in the total angular momentum squared due to Euler and Adams- 
Bashforth second-order integration are, respectively, 


and 


AN n = h 2 u n 2 N n sin 2 x n 


(87) 


AN n = 


— T-i. 

-hL n A L 


n-1 



(88) 


where Xn * s the an S* e between the c«J n and L n vectors. Again, the change due to 
Euler integration is positive definite and the Adams-Bashforth integration result is not 
positive definite. 

Near cylindrical symmetry: Assume that the X-axis to be the axis of near sym- 
metry, that is 


*33 = *22 + 51 ( 89 ) 

By using definition (89), equations (86) become, if first-order terms in 61 are retained: 

Li = - -iL l 2 l 3 
* 22 2 


L2 

l 3 



KLiL 2 


L 3 Ll 




(90) 


where 
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The near cylindrical symmetric approximation is used when 


for which equations (90) degenerate 
tions (90) become, in matrix form, 


K l 

1 61 


I <<c l 22 


1 1 V 

to the equations 


of an oscillator since 


Ll « 0 . 


with 


^2 


0 

KLi 


l 2 

^3 


-KLi 

0 


L 3 


KLj * K3np 0 




(91) 


Equa- 


(92) 


(93) 


where p Q is a constant roll rate. 

To analyze equations (92), the range of the quantity Kin needs to be determined. 
To accomplish this end, first observe that the inertias do not uniquely define the body 
configuration; that is, bodies of uniform material but with many different shapes can all 
have the same inertia matrix. Hence, perturbations from the set of all cylinders of uni- 
form material are (dynamically) equivalent to the set of all nearly cylindrically symmet- 
ric bodies. Therefore, no generality is lost by considering perturbations of the cylindri- 
cal body of uniform material shown in figure 9; the inertias are (for 61 = 0): 

Ill=|mR 2 ^ 


I 2 2=^m(3R 2 + J 2 )> 


J 33 = l 22 


(94) 


Where R and l are as indicated in figure 9. Substituting equations (94) into equa 
tion (93) yields 
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< 9! 

3 + r^ 

where r = l/R is the ratio of length to radius and always positive. Looking at the 
derivative of K with respect to r 

I,, dK , < 0 

dr 3 + r 9 

reveals Kill to be a monotonic decreasing function of r (r > 0). The minimum and 
maximum values of IllK are at the maximum and minimum values of r. Looking at 
the limits of the rotating disk (r — 0) and rotating rod (r — °°) establishes the limits of 
equation (95). 


lim KLi = p 


lim KLj = -p 

r-*°o 


The general character of InK is shown in figure 10. Consider a slight perturba- 
tion from cylindrical symmetry 

J 33 = l 22 + 61 

For the near cylindrical symmetric approximation to be valid 


122 K 


but from equations (91), 


J 22 " J 11 


Hence, as shown by this relation, small perturbations from cylindrical symmetry for very 
long slender and/or short flat bodies do not greatly remove the motion from that of cylin- 
drical symmetry since relation (97) is valid. However, as 


38 



the motion approaches a region where nonlinear motion resumes since the coefficients 

~j6l[ 


in equations (90) are then on the same order of magnitude that is, K 


Hence, 


l 22 


figure 10 can be divided into two regions. Region I is where relation (97) is satisfied and 
the motion is nearly linear. Region n is where relation (97) is not satisfied and the 
motion is not described by the approximate equations (92). In region II is the class of 
bodies with near spherical symmetry. 

For near cylindrical symmetry, the angular momentum rate equations are 


(98) 


i-2 


0 -k 


L 2 



k 0 


L3 


where k = LjK. 


Consider the integration of equations (98) by using first the Euler integration and 
then the Adams-Bashforth second-order integration formula. For Euler integration, 




L2 1 


0 

-k 


L2 

l 3_ 

n+1 

L3 

= h 
n 

k 

0 


l 3_ 


Taking the z-transform and collecting yields 


(z - 1) - h 

which has a nontrivial solution if 

(z - 1) - h 

The characteristic equation of equation (100) is 

(z - l) 2 + h 2 k 2 = 0 



0 -k 


Li(z) 

k 0 


L2( z ) 


(99) 


( 100 ) 


( 101 ) 


( 102 ) 
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Equation (102) can be factored into two roots by inspection; the roots are 


z - 1 ± ihk 


(103) 


from which the solution to equations (99) by using approximation (35) is 


L3 


hie 


= e 




-ik.i¥t 


+ c 2 e 


(104) 


where cj and C2 are constant vectors. By using equations (104), the solution to 
equation (87) is 


■vt t ^ hk^t 

N = L l,o +e 


L 2,o +L 3,o 


(105) 


Note that for the case of cylindrical symmetry, the rate equations are linear (see 
eqs. (98)). For this case an expression can be found for the error when Adams-Bashforth 
second-order integration is used. Consider 




l 3 


n+2 


L2 


l 3 


n+1 


0 -k 


k 0 




l 3 


n+1 L 


0 -k 


k 0 


l 2 


l 3 


n 


(106) 


Taking the z-transform and collecting terms yields 



0 -k 


L 2 ( z ) 

z(z - 1) - |(3z - 1) 

k 0 

L Jj 

/ 

_L 3 (z) 


= 0 


(107) 


for which the characteristic equation is 
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z 2 (z - l) 2 +ll^(3z - l) 2 = 0 
4 


(108) 


which can be factored as 


z(z - 1) = ±i ^(3z - 1) 


The roots of this equation are known to be 

Z 1 2 s 1 - 2c 2 - 2c^ ± i^2c + 2c^) 

z 3 4 ~ 2c ^ + 2c ^ ± i( c “ 2c ^) 


(109) 


(HO) 


(111) 


hk 


where c = — . (Compare eqs. (57) and (58).) The roots given by equation (111) are 

heavily damped and need not be considered. The solution obtained by using approxima- 
tion (35) is 


L3 


,3,4 

4 h_J^ t 

, 16 


f 1 ( k h2k3 ) t _i ( k + 12 h2k ^ 

j C G + C2^ 


c^e 


J 


( 112 ) 


where cj and C 2 are the same as in equations (104). The solution to equation (88) is 
given as 


8 h°k* t 

M _ T 2 L 2 T 2 ] 16 

N “ L l,o + ( l 2,o +l 3,o J e 


(113) 


Near spherical symmetry: Consider region II of figure 10 where nonlinear effects 
cannot be neglected (that is, relation (97) is no longer valid). In the nonlinear region K 
is no longer large compared with Sl/ 122 2 ; this condition implies that 


r 22 " *11 

hi 


J 33 " J 22 
122 


(114) 
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A solution to the general problem has not been found in region II. However, an approx- 
imate solution can be found for the Euler integration formula if, in addition to rela- 
tion (114), the following relations also hold: 

I I "N 
6I t 


111 

61 


« 1 


« 1 


(115) 


122 


where 

61 ,= I 22" I 11 ( 116 ) 

61 = I 33 - I 22 ( 117 ) 

Relations (115) then define near spherical symmetry and equation (12) can be rewritten as 


— — 
p 


Ll 


q 

_ 1 
!22 

l 2 

+ 

r 


L 3 







61 ' 


0 


0 


122 

0 0 
0 0 


61 


l 22 


Ll 


L2 


L2 


+ O 


er 

V 22 '’ 


(118) 


Substituting equations (118) into equation (87) and neglecting terms in 61 and 61 
results in 

h2vr 2 

h N n 2 
AN n * — sin z Xn 

I22 2 

Dividing equation (119) by N n ^ and letting 

h 2 . 2 

«n = o Sin X n 

122 


(119) 


( 120 ) 
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yields 


AN 


n 


- K n 


N. 


n 


( 121 ) 


For AN n « N n , equation (121) can be approximated by a perfect difference, namely, 


M aN n - AN n 
W ' N„N n+1 Nn 2 


( 122 ) 


The approximate solution to equation (88) for nearly spherically symmetric vehicles when 
approximation (122) is valid is 


N r 


N r 


1 - N 0 2 “i 
j=o 


which for /cj 


nearly constant or constant 


N(t) * 



and the resulting error in the total angular momentum squared is 


e(t) * 


No 


Npttt 
h - N 0 fct 


(123) 


(124) 


(125) 


Thrusted case: To study the on-off thrust controlled vehicle, assume that M is 
defined only on the integration intervals. This assumption neglects the small delays in 
applying thrust. The effects of this delay can be evaluated by analysis or by comparing 
solutions at two different interval sizes. The latter method was chosen since the fine 
detail of the limit cycle motion was not of great interest in this problem. Under this 
assumption, the M of equation (11) can be integrated exactly by using Euler integration. 

By using Euler integration, the per-step change in the total angular momentum 
squared is 

AN n = hN n + h 2 (u) n 2 N n sin 2 x n + 2co n N n 1//2 M n sin x n cos y n + M n 2 ) (126) 
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where M is the magnitude of M and y is the angle between the L and M. Nom- 
inally, the control moment M is applied to w. Since L and u> are orthogonal, y 
is always near 7t/2. Hence, 

AN n - hN n + h 2 (w n 2 N n sin 2 x n + M n 2 ) (127) 


From equation (127), the per-step deviation from Euler integration is 


Ae n 



(128) 


The importance of this result is that the per-step error given by equation (127) is inde- 
pendent of the frequency associated with M; only the magnitude M contributes to the 
solution error. 

It would seem that an impasse has been met since on-off control thrust dictates the 
use of Euler integration. However, equation (125) indicates rather poor behavior of the 
nonthrusted portion of the moment equations for this scheme. Equation (11) can be par- 
titioned so that each term can be integrated in a manner most likely to yield a quality 
solution. To do this, consider the integral of equation (11). 


L = 


I( n ■ L >“ + f 


M dt 


(129) 


Define the partitioning as 


Lj = . h)dt 

L 2 = Jm dt 

L = Li + L2 


(130) 


The two terms, Lj and L 2 , can be integrated by using the Adams-Bashforth second- 
order and Euler integration schemes, respectively, and a high quality solution is obtained 
with little added complication. 

Discussion of exper im ent .- The angular momentum rate equations were first 
analyzed for the nonthrusted case ^that is, M = o) for two integration schemes. The 
analysis was done for two classes of vehicles which were nearly cylindrically symmetric 
and nearly spherically symmetric. The forms of the angular momentum rate equations 
are the same as those of the direction cosine rate equations with the exception of the 
nonlinearity defined by uJ as a function of the L. For the near cylindrical case, the 
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equations are nearly linear and the solution for Euler and Adams-Bashforth second-order 
integration can be approximated. 

The angular momentum rate equations were solved for initial conditions of 

D = a = r _ 1 radians 
^3 sec 

111 = 500.0 kg-m^ 

122 = 878.9 kg-m^ 

I 33 = 881.9 kg-m^ 

h = — sec 
32 

with Euler and Adams-Bashforth second-order integration schemes. For these values of 
inertias, the motion is described by the near cylindrical approximation. The predicted 
and computed motion for Euler integration is shown in figure 11 and the agreement is 
good. The corresponding predicted and computed errors for Adams-Bashforth second- 
order integration is shown in figure 12. Two computed errors are shown in figure 12; 
they are for Nj-, = 27 bits and = 48 bits. The roundoff error is shown to be domi- 
nant on the 27 significant bit machine. The solution using extra precision is closely 
approximated by the approximate analytic solution of the Adams-Bashforth second-order 
integration. It is not clear at this time whether the roundoff error is caused from the 
induced error or the inherent error. The Adams-Bashforth scheme is sufficient for 
solving the nonthrusted portion of the angular momentum equations, the accuracy 
exceeding that of Euler integration by several orders of magnitude, even for = 27 bits. 

The use of Euler integration for a vehicle of near spherical symmetry is, in gen- 
eral, inadequate. The behavior of the total angular momentum squared was found not to 
be a positive definite form under Adams-Bashforth second-order integration. In general, 
one can expect the error due to Adams-Bashforth integration to be of several orders of 
magnitude less than that due to Euler integration. 

It was assumed that the small time delays in applying thrust were negligible. This 
assumption clearly depends on the control system parameters. For the Gemini problem, 
the effects of these time delays were determined by comparing two solutions of different 
interval sizes which were h = 31.25 milliseconds and h = 15.625 milliseconds. The 
control system was placed in attitude hold and three stick deflections were commanded 
as follows: 
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(1) Roll right at one -half maximum rate for 0 = t = 3 seconds 

(2) Pitch up at one-half maximum rate for 10 = t = 13 seconds 

(3) Yaw right at one-half maximum rate for 20 = t = 23 seconds 

The results of these attitude -control cases are shown in figures 13 to 15. As seen from 
the figures, the gross behavior is, for all purposes, identical. There are small apparent 
discrepancies in the limit cycle and end point motion. Even so, these solutions show 
good agreement. 

When L and To are parallel, the integration of the nonthrusted portion is 
integration-scheme independent. The L and To are parallel whenever the To ^hence, 
the l) lies along a principal axis. Most vehicles are constructed so that the body axis 
and principal axis are nearly alined. The typical pilot uses the single-axis technique 
(that is, a sequence of maneuvers about a single body axis at a time) for which To and 
L are nearly parallel. The primary contribution to the angular momentum is thrust 
which can be integrated exactly. 

When either the L and To are nearly parallel and/or the magnitude of To is 
small, the roundoff error is the dominant error. Clearly, the most likely area of diffi- 
culty is in the limit cycle motion. The results indicate no serious instability in the limit 
cycle motion; thus, the only remaining question is the fuel usage in the limit cycle. The 
change in fuel usage will not be greatly affected for short-period flights; however, for full 
mission simulators where the motion is in the limit cycle for several hours, considera- 
tion should be given to the fuel cost due to computation errors. 

In simulating the angular momentum rate equations on analog computers, the mini- 
mum pulse width of the attitude control moments is determined by the switching time of 
the computers logic components unless otherwise provided. An analogous situation 
occurs in digital simulation where the minimum pulse width (Euler integration) is the 
integration interval size unless it is provided by other means. In digital simulation prob- 
lems where rather large integration intervals are used, a minimum pulse width must be 
provided. Large erroneous expenditures of fuel due to erratic limit cycle motion can be 
caused by oversized minimum pulse widths which are determined by the integration 
interval. Although the simulation of minimum pulse width on general purpose analog 
computers is somewhat difficult, this provision is an inherent consequence of digital 
simulation and requires only a small amount of simple logic. 

It is often desired to compute the mass change of the simulated vehicle. The mass 
rate is usually given as proportional to thrust. Hence, Euler integration gives an exact 
representation of the mass since mass is a linear function of time or a constant. 
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Translational Dynamics 

Analytical development .- The mathematical model representing the trajectory 
motions is extremely nonlinear and therefore does not lend itself very readily to the type 
of analytical investigation made in the other sections. However, since the trajectory of 
the two vehicles referenced to inertial space is very smooth, some comments concerning 
truncation compared with roundoff effects can be made. The frequencies associated with 
the trajectories are on the order of 1.2 x 10"3 radian/second. When this frequency is 
coupled with the maximum interval size, h ~ 0.05 second, the graph in figure 2 provides 
an estimate of the local truncation and roundoff errors (for hco = 6 x 10~5 radian): 

e tr a 6 x 10“ 10 
e r « 6.2 x 10 -5 

The indication given by these results is that roundoff error constitutes most of the error 
in the solution of the trajectory motion. In order to verify the prime source of error and 
to evaluate the effects of the propagated truncation and roundoff errors, computer solu- 
tions with widely varying interval sizes can be compared with a double-precision fourth- 
order Runge-Kutta solution (or any other reliable solution). 

Discu ssion of experiment .- The global errors, using double-precision fourth-order 
Runge-Kutta as a standard, in the trajectory variables (eqs. (Al) to (A16)) are shown in 
figure 16 as functions of the integration interval size at t = 5004 seconds. It is noted in 
figure 16 that the minimum global errors occur around an interval size of h ~ 0.1 second. 
However, because of the pilot input and the need for a smooth visual display, an upper 
bound on the interval size must be placed at 0.05 second (that is, h % 0.05 sec). The 
logical choice, if constrained to powers of two to minimize roundoff effects, is 
h = 2 - ^ second. Also, note that there is little difference in the global errors of Euler and 
Adams-Bashforth second-order integration; thus, truncation error constitutes a minor 
part of the total error and hence, contributes to the confirmation of the validity of error 
dominance. In figure 17 is the global error in two coordinates as a function of time for 
h = 2~5 second. Again, it is apparent that roundoff error is by far the largest contrib- 
uting error. 

Further evidence verifying roundoff error (in particular, the induced error) as the 
major contributor to the global error is seen from the results in table H. There is excel- 
lent agreement between the standard solution and the partial double-precision (for control 
of induced error) second-order Adams-Bashforth solution. When compared with the 
single precision second-order Adams-Bashforth solution, the severity of the roundoff 
effects on a = 27 bit class of machine is clearly indicated. 
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TWO EXAMPLE PROBLEMS 

In the past, the simulation engineer has relied heavily on an intuition developed 
from years of experience with analog (electronic differential analyzer) computation. The 
advent of digital simulation requires the building of a reliable intuition of the numerical 
processes of digital computation before digital simulation can become highly productive. 

It is for this reason that two example problems are considered. 

The first problem is a Gemini-Agena rendezvous and station-keeping simulation. 
(See ref. 11.) This problem fits the class of problems studied in the analysis. Hence, 
the integration schemes discussed are directly applicable. 

The second example is a modification of the first, the modification being the addi- 
tion of an elastic tether (ref. 12) coupling the two vehicles' motion. The two vehicles are 
now somewhat like a single structure with internal degrees of freedom. Considered as a 
single vehicle, the system has a structural frequency that depends on the tether proper- 
ties. The results of the analysis are no longer applicable in its entirety. However, 
these results are used as a starting point. The methods used to improve the digital pro- 
gram are discussed later. 

Hopefully, several things will be accomplished in these two examples. FORTRAN 
is adequate in its present form for the IBM 7094 -H, especially with regard to the amount 
of core and processing time for large simulation problems. The storage and processing 
time depend upon how well the FORTRAN compiler has been adapted to the particular 
machine being used. The integration schemes discussed in the analysis were not intended 
for blanket approval. The second example problem then gives the next logical step for 
arriving at a working simulation program. The tether program is a direct modification 
of the first. Hence, the effort involved to modify a large existing digital simulation pro- 
gram is better understood. 


Gemini-Agena Problem 

Program developme nt.- The physical problem being simulated was that of two 
vehicles in near-earth orbit. One of the vehicles was constrained to a planar orbit and 
was considered to be passive. Tumbling or any special orientation could be simulated 
with proper initial conditions on the inertial moment equations. No external thrust was 
provided. The other vehicle had a full six degrees of freedom. There were external 
thrust capabilities provided on this vehicle. The external thrust capabilities can be 
divided into two categories. The first category includes translational movement jets 
which were located approximately along each principal axes. These jets were used in a 
strict on-off fashion and were activated by the pilot's stick deflection. External thrust 
for the second-category attitude movements can be further subdivided into several modes 
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of operation. These jets could be activated by a system of switches and stick movements 
in either rate command with attitude hold, acceleration command, or pulse modes. Each 
of these modes was determined by using standard methods. A block diagram depicting 
the flow of the equations (eqs. (A4) to (A106)) is presented in figure 18. The control sys- 
tem was not listed in appendix A because of the difficulty of expressing the control sys- 
tem in equation form. 

An exception has been made to the original Gemini-Agena problem. The "cage” 
mode of the Gemini control system has been deleted because of the present limited num- 
ber of digital-analog converters but this deletion does not jeopardize the results of this 
study. 

D iscussion of experiment .- The Gemini-Agena station keeping and docking problem, 
the equations of which are discussed in appendix A, were programed for the IBM 7094 -EL 
in FORTRAN. This program was used with the Langley visual rendezvous simulator 
which is shown schematically in figure 3 of reference 12. Some of the results of this 
program are discussed in references 11 and 12. The integration schemes used in this 
program are given in the following table: 


Portion of program 

Integration scheme used 

Euler parameters 

Adams-Bashforth second order 

Control system 

Adams-Bashforth second order 

Angular momentum: 
Dynamic part 
Thrust part 

Adams-Bashforth second order 
Euler 

Trajectory: 
Dynamic part 
Thrust part 

Adams-Bashforth second order 
Euler 

Mass 

Euler 


Subroutines such as sine, cosine, and square root provided with this language were used. 
This Gemini-Agena program was used in conjunction with a machine language program 
which provided mode control and communications with a real-time clock and simulation 
hardware. The integration interval used was h = 31.25 msec. The two programs com- 
bined required about 6000 locations in core and less than 4 msec for all calculations and 
input/output functions. This program leaves approximately 85 percent of core unused 
and the central processor free approximately 87 percent of the time. 

The amount of core and central processing time required for this problem clearly 
indicates that more than one such problem could easily be processed simultaneously 
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(multiprocessing) with sufficient conversion equipment and a control program to provide 
independent mode control. (See ref. 3 for other discussions of multiprocessing.) This 
operation is limited in that a single program must contain all such problems to be com- 
piled as a group. On a more sophisticated machine, it would be desirable to compile, 
process, and terminate a real-time simulation job while processing other real-time pro- 
grams and/or nonreal-time programs whenever sufficient core or central processing 
time is available (multiprograming). 

Gemini -Agena Tether Problem 

Program development .- The system consfdered is essentially a modification of the 
Gemini-Agena rendezvous and docking problem already considered. The equations of 
motion are referenced to the center of mass to retain a similar form for the translational 
equations except, of course, for the elastic coupling. A more complete discussion of 
this problem is in reference 12. 

The longitudinal translational vibration frequency was estimated to be 0. 5 Hertz for 
the worst case. The expected relative local truncation and roundoff error for this mode 
is (o>h ~ 9.6 x 10 - 2 radian) 

e tr * 1.5 x 10“ 4 
* 3.8 x 10“® 

for h = 31.25 msec where a second-order integration technique is assumed. The new 
equations of motion were programed and Adams-Bashforth second-order integration 
schemes were used for both the translational and rotational degrees of freedom. The 
program was initially checked by using this method of integration and the answers were 
compared with a fourth-order Runge-Kutta solution. The observed computation errors 
(Error = | Runge-Kutta — Adams-Bashforth I/maximum value) were on the order of or 
less than 10 percent. 

At this point a simple analysis was made. The initial conditions can be chosen so 
that the longitudinal translational equation becomes 

.. 9 

X = -<jO“X 

where co is related to the masses of the vehicles and spring constant of the tether. The 
roots using Adams-Bashforth second-order integration are found to be 

Zj g ~ 1 - 2c^ - 2c 4 ± i2c(l + c2) 
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z g 4 ~ 2c 2 + 2c^ ± ic(l - 2c^) 


where c = hco/2. 


yields 


Transforming these roots to the s-plane (s log^ z) and using approximation (35) 

1 h e 


Sj g * h 3 o> 4 ± iw fl + ^ h 2 co2j 


with the real parts of Sg ^ large and negative. 


The next scheme considered was to use Adams-Bashforth second-order scheme to 
obtain x and the Adams-Moulton second-order scheme to obtain x: 


*i+2 = *i+l + \ ( 3x i+l - x i) = *i+l - ( 3x i + l - x ij 


and 


x i+2 = x i+l + §(*i+2 + x i+l 


which can be written in matrix form as 


r~ 

' 


— — ■< 



E(E - 1) 

h|!(3E - 1) 


x i 


0 

- | E(E + 1) 

E(E - 1) 


1 

X 

1 ^ 


0 


If the z -transform (E — z) is taken, the characteristic equation is 


z(z - 1)^ + c2(z + l)(3z - 1) 


= 0 


Hence, one root is zero and the factor in brackets can be written as 


(z - c^) 


z^ - 2(1 - 2 c2)z + 1 


+ 4zc^ = 0 


If the system is assumed to be nearly stable ^that is, | z | % l) , the last term is 


4zc 4 ~ 10 


-6 
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Si 9 ~ - — ± ico f 1 

where Sg is large, real, and negative and s^ = - eo . Hence, the roots are by a factor 
of four closer to the real roots in damping and the frequency error has improved by 
almost a factor of two. 

Discussion of experiment .- The Gemini -Agena program, discussed in the previous 
section, was modified to accommodate the tether and the additional degree of freedom. 

The integration schemes were not changed in the initial stages of checkout. The observed 
computation errors were on the order of or less than 10 percent when compared with an 
all-fourth-order Runge-Kutta solution. The Adams-Bashforth scheme was then modified 
by applying the Adams -Moulton second-order corrector once to all translational and 
rotational variables. The solution obtained from this integration scheme was within one- 
half percent of the fourth-order Runge-Kutta solution and the addition to the computation 
time is virtually undetectable. The final integration schemes used are summarized as 
follows: 

Integration scheme used 

Adams-Moulton predictor-corrector* 
Adams-Bashforth second order 

Adams-Moulton predictor-corrector* 

Euler 

Adams-Moulton predictor-corrector* 

Euler 

Euler 

*Second order with corrector applied once. 


Portion of program 

Euler parameters 

Control system 

Angular momentum: 
Dynamic part 
Thrust part 

Trajectory: 

Dynamic part 
Thrust part 

Mass 
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The amount of core and central processor time required was virtually no different than 
that required for the Gemini-Agena problem. The modification and ensuing checkout 
process was accomplished in less than one man-day (8 hours). This time does not 
include the Runge-Kutta check case that was already available. 

RESULTS AND DISCUSSION 

The basic equations considered in the analysis are common to most flight simula- 
tion problems. In particular coordinate transformations, the attitude control system and 
angular momentum equations commonly reoccur. The greatest variations are in the 
addition of aerodynamics, structural effects, and the description of the trajectory which 
could cause large deviations from the conclusions of this study. For example, the addi- 
tion of nonlinear aerodynamic moments could greatly alter the basic stability of the 
angular momentum equations as discussed in the analysis. In some special applications, 
the trajectory is computed from the body-axis force equations which are of the form of 
the direction cosine rate equations, the constant of motion being the total velocity squared 
for no applied forces. The addition of body bending modes and the flutter of appendages 
generally raise the frequencies of the problem. The effects of truncation error must 
then be reexamined and a higher order integration scheme may need to be applied. 

The trajectory calculations of many simulations, in general, exhibit a low-frequency 
component, where the frequency of this component is several orders of magnitude smaller 
than the frequencies associated with the remainder of the problem. Hence, roundoff error 
will generally be a concern. The analysis of roundoff and truncation error should be 
helpful in determining the need for concern. Hence, if a particular problem has the local 
relative roundoff and truncation error of the same order of magnitude, great care should 
be applied to determine the propagated roundoff error. If roundoff error is important, 
partial double precision should be applied to as few of the sensitive variables as possible 
to obtain a reasonable solution. 

Roundoff error has been shown to be a prime consideration on a computer with a 
27-bit fractional part. Langley Research Center now has available a complex of Control 
Data series 6000 computers. The fractional portion of the floating-point word is 48 bits 
in this computer series; hence, the roundoff error is reduced by 2~21 as compared with 
the error for the N^ = 27 bit machine used in this study. 

Because of the roundoff error, the choice of integration interval size is greatly 
limited for manned simulations of this class of problems on an Nfo = 27 bit computer 
^that is, 31.25 msec •g.h = 50 msec where the lower limit is chosen to control roundoff 
error). Greater flexibility is afforded on the 6000 series machines (namely, 

2~21 x 31.25 msec ~ 10 - 5 msec ^h = 50 msec) without roundoff error becoming a serious 
problem. 
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APPLICATIONS 


Langley Research Center now has three Control Data series 6000 computers with 
multiprograming capability in a mixed real-time and/or nonreal-time environment. 

The prime computer for real-time jobs is a 6600 computer with a 131K)jq memory. 

Most of the problems solved digitally to this date are problems which are clearly too 
large or marginally large for the analog-TRICE computer complex (TRICE, a digital dif- 
ferential analyzer) at the Langley Research Center. There are several real-time simu- 
lation consoles which provide the simulation engineer program control, input, and output 
(to change problem parameters). The language used with this computer complex is the 
FORTRAN language with a relatively early compiler; the compiler is still relatively 
inefficient with regard to central processor time. A brief description of several flight 
simulation problems with corresponding vital statistics that have been completed or are 
currently in progress follows: 

The lunar-orbit and landing-approach simulations were formulated for real-time 
digital operation. This problem consists of three separate vehicles each employing six 
degrees of freedom and relative geometry between vehicles. All integrations used 
second-order Adams-Bashforth and Euler integrations as previously discussed. The 
results were largely compared with analytic results and showed good agreement. The 
interval size was 1/64 second or 15.625 milliseconds (this size could be increased) which 
was chosen to accommodate the minimum impulse of the attitude-control system. The 
required central processor time was approximately 7 milliseconds including all overhead 
The memory requirement was approximately 37K)^q. (37K)iq denotes 37 000 in 

decimals.) 

An Apollo-LEM rendezvous program was obtained by modifying the lunar-orbit and 
landing-approach program. The integration schemes and interval size were the same. 
The central processor time was approximately 8.5 milliseconds with 24K)jq core loca- 
tions required. 

The lunar-orbit and landing-approach program was again modified for an Apollo- 
abort study. Among the features retained in this modification is staging capability. The 
central processor time was approximately 7 milliseconds and less than 24 K)iq core loca- 
tions were required. 

The original Gemini-Agena program, the first example problem for the IBM 7094-n 
was converted for the Control Data 6600 computer. The interval size was 2-5 second 
or 31.25 milliseconds when second-order Adams-Bashforth and Euler integrations were 
used. The central processor time was approximately 2 milliseconds with 6K)io core 
locations required. 
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A simulation program of F-105 and F-86 fighter airplanes has been developed to 
study the feasibility of employing research simulators for evaluating the tactical effec- 
tiveness of fighter aircraft. In the program, each airplane had six degrees of freedom, 
relative geometry equations, and simplified line-of-sight scoring equations. Adams- 
Bashforth second-order integration was used exclusively with an interval size of 2 - ^ sec- 
ond or 31.25 milliseconds. The error was less than 0.5 percent when compared with 
fourth-order Runge-Kutta check cases. The total central processor time including 
input/output and the real-time monitor is approximately 17 milliseconds and requires 
less than 4 OK) jq core locations. (Note that this program includes two complete airplane 
simulations.) 

The HL-10 lifting body was simulated in full six degrees of freedom using second- 
order Adams-Bashforth integration with an interval size of 2“5 second or 31.25 milli- 
seconds. The errors were again under 0.5 percent. The total central processor time 
was 10.21 milliseconds with 16.5 K)iq core locations required. 

The static test program is a FORTRAN- coded simulation of an orbiting laboratory 
of an Apollo telescope mount (ATM) configuration. Six elastic-body modes are incorpo- 
rated directly in the model, the other modes being accounted for in a quasi-static model. 
The program is being used to evaluate the use of a control-moment-gyroscope system to 
generate control torques for attitude stabilization of the telescope mount. This program 
runs in a closed loop with three control-moment gyroscope prototypes mounted on a 
moving-base simulator. The loop is closed through 14 bits plus the sign, analog-digital- 
conversion equipment which transmits the measured torque output of the control-moment- 
gyroscope system. There are provisions for external torque input from taped frequency- 
modulated analog or discrete crew motion disturbances. The fourth-order Runge-Kutta 
integration scheme is used with h = 31.25 milliseconds which is sufficient to keep errors 
below 0.2 percent when compared with a nonreal-time solution. The central processor 
time is about 19 msec with approximately 16 K)iq core locations required. 

In addition to these large simulation programs, many small simulation problems 
have been studied. Many advantages have become apparent from these small jobs. These 
jobs can be used to fill central processor time slots and core positions not used by the 
large simulation jobs (which is otherwise filled with background nonreal-time jobs). Many 
research engineers are familiar with the FORTRAN language; hence, the research engi- 
neer can prepare his own simulation program without relying on highly specialized per- 
sonnel. It has been demonstrated that inexperienced personnel can operate even large 
digital simulation problems with little supervision. 
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CONCLUDING REMARKS 


The analysis indicates that the lower order integration schemes, second-order 
Adams-Bashforth and Euler, are generally adequate for the class of nonaerodynamic 
rigid-body problems. The accuracy of these techniques indicates that these integration 
schemes are reasonable starting points for similar flight simulations. In the event of 
failure, a simple analysis, similar to those contained in this paper, can often be performed 
to indicate the next step in determining a working simulation program. The Adams- 
Moulton second-order corrector is a reasonable second step. A relatively efficient com- 
puter program in either case is indicated by the computer results of this study. 

The Green's function technique employed to solve the gyroscope equations is exact 
when the forcing function is linear over the integration interval. Higher order terms in 
the forcing function can be carried in this technique with effectively no added computer 
time. This method should be extended to include such effects as nonwindup limiting 
(piecewise linearity). This technique can also be extended to other linear and piecewise 
linear systems. For example, consideration should be given to the computation of struc- 
tural effects where frequencies are often high. 

Assuming the control moment to be defined only on the ends of the integration inter- 
vals and hence neglecting small time delays in applying thrust was shown to be valid for 
the Gemini-Agena example problems. This procedure was also valid for other vehicles. 
Clearly, this approximation is dependent on the control system parameters and further 
analytical work is needed. However, a simple check by changing the integration interval 
size should lend confidence in the solution of any particular problem. 

Roundoff error was shown to be marginal on a floating-point machine with a 27-bit 
fractional part with an integration interval of 31.25 milliseconds. Partial double pre- 
cision was shown to be an effective method for controlling the induced error. It is 
recommended that double precision be used with care because of the additional central 
processor time. On machines with a fractional part less than 27 bits, partial double pre- 
cision is likely to be unavoidable. In any case, efforts should be made to control the 
inherent error and error dominance should be used on short word machines (less than 
27-bit fractional part) to locate areas where the induced error is likely to cause difficulty. 

Low-order multistep integration techniques were applied to two large simulation 
problems and solved on an IBM 7094-H (a large last-generation digital computer). The 
results of these experiments show these problems to be small compared with the capacity 
of this computer. Hence, multiprocessing (concurrent processing of more than one prob- 
lem but in a dependent manner since all problems must be compiled as a single program) 
on the last-generation computer is a reasonable approach to reduce cost. Otherwise, a 
smaller slower machine (hence, less expensive) is sufficient for many simulation problems. 
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Recent computers with hardware and software oriented toward multiprograming in 
a real-time environment (concurrent processing and/or compiling more than one problem 
independently) are clearly indicated as desirable for a general-purpose simulation labo- 
ratory. Langley Research Center now has a complex of Control Data series 6000 com- 
puters with a multiprograming system. This mode of operation has proven to be very 
successful in that one or more large real-time simulation programs can be processed or 
compiled, with other smaller real-time or nonreal-time jobs giving better utilization of 
the computer's capabilities. 

The use of FORTRAN as a simulation language is shown to be adequate, especially 
in view of the core storage and central processor time required for the several simula- 
tion problems discussed, most of which involved two or more vehicles. However, the 
efficiency of the FORTRAN language depends on how well the compiler has been adapted 
to the computer used. The efficiency of the object code of the 6600 computer used for 
this computation is not yet commensurate with the capabilities of this machine. 

This investigation has given primary consideration to the simulation of nonaerody- 
namic vehicles. The conservative attitude about possible instability when aerodynamic 
control is added to the simulation is largely unwarranted. Two rather unconventional 
aircraft simulations, the F-105 airplane and the HL-10 lifting body, indicate that the inte- 
gration schemes are applicable to a wide class of aerodynamic problems. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 21, 1968, 

125-23-03-02-23. 
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APPENDIX A 


GEMINI- AGENA EQUATIONS OF MOTION 

By Roland L. Bowles 
Langley Research Center 

The purpose of this appendix is to present in concise format, with no derivation, 
the basic equations used in this study. The following equations have been developed to 
simulate the Gemini- Agena configuration in eleven degrees of freedom. These equations 
are particularly well suited for simulation of terminal rendezvous and such close-in 
operations as station-keeping and fly-around missions. Foremost in the selection of 
this formulation was that it meet mission and simulation hardware (visual docking 
simulator) requirements and overcome computing difficulties with this class of piloted 
simulations. 


Equations of Motion 

The force equations defining relative motion of the observer vehicle (three degrees 
of freedom) are referred to a rotating local vertical axis system centered in the target 
vehicle (two degrees of freedom). Polar coordinates (r s ,0 s ) are employed to locate the 
center of gravity of the nonthrusted target vehicle with respect to inertial space. There 
is essentially no restriction imposed on the eccentricity of the target orbits. Basic coor- 
dinate systems and relative geometry of the two vehicles are shown in figure 1. Both 
vehicles are assumed to be moving over a homogeneous nonrotating spherical earth. 

In order to improve absolute computational accuracy and increase effective com- 
puter resolution, perturbated equations were used to generate the target orbit. For tar- 
get orbits of moderate eccentricity, there are obvious computational advantages to be 
gained by transforming the dependent variables to represent deviations from a nominal 
circular orbit. In this manner the importance of computer errors for orbital or near- 
orbital simulations can be minimized. The necessary dependent variable transforma- 
tions are 



v 0 = l/g 0 r o 


(Al) 


(A2) 
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where r Q is the radius of the reference orbit, g Q is gravity acceleration at r s = r Q , 
and V 0 is the characteristic velocity of a vehicle in a circular orbit at radial distance 
r Q . In equation (A2), V is defined as the tangential velocity at any time (that is, 

V = -r s 0 s ). For computer scaling purposes, it is advantageous to choose r Q equal to 
the mean radius of the highest apogee and lowest perigee trajectories which must be 
computed. Hence, 


r o = r e 




(A3) 


where r e is the earth radius and h a and hp are the apogee and perigee altitudes, 
respectively. Thus, the nondimensional perturbation variable p will have a symmetric 
excursion about p = 0, over one orbital period. The degree of symmetry for 6 depends 
strongly on the eccentricity of the orbit to be generated. By using the definitions of 5 
and p in addition to the well-known Keplarian equations of motion for two-dimensional 
orbits and by making use of the fact that the orbital angular momentum is a constant of 
the motion, the resulting perturbation equations can be written as 


1 Kq 2 + 2K q-p 

C0 0 2 (1 + p) 3 


(A4) 


8 


s ~ 


-OJ 


o 


1 +Kq 

_(i + p) 2 _ 


where cu 0 2 = and Kq 
Kq = 5(0) [l + p(0j] + p(0)). 


is specified from initial conditions (that is, 

The coordinates (r s ,# s ) are computed as follows: 


r s = r 0 (l + P) 


0s 



9 S dr 


(A5) 


(A6) 

(A7) 


Equations (A4) and (A5) are characterized by divisions which have a denominator of the 
form (1 + p) n where n is an integer and typically p « 1. 

The observer vehicle is located relative to the target local vertical by the three 
components of R that is, (r = x,y,z). Since |r| « r s for missions considered in 
this simulation, the solution can be considerably enhanced by expanding the gravity 
gradient which appears in the equation for R and by retaining terms to the desired 
order of accuracy. Relative equations which are correct to the second order in the 
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gravity gradient were utilized in this study. The scalar equations used to define relative 
motion of the observer vehicle are 


x + 20 s z + 0 s z 



3a>o^ z 
r o (1 + p)3 


m 


(A8) 


.. . . .. • n 2w n ^Z 

z - 26 s x - 6 s x - z — - , 

S S (1 + p)3 2 r Q (1 + p) 4 


+ 3 u o* 


(x% + y2 




(A9) 


y + 


Wq2 3o> 0 2 yz 

(1 + p) 3 y r o (1+p) 4 



(A10) 


where T-jj, Ty, and Tg, are components of observer vehicle control forces resolved 
in the target local vertical axis system, and m designates the constant mass of the 
observer vehicle. Acceleration coupling terms 6 S and p resulting from elliptical 
target motion are retained. The variables 6 S , 6 S , p, and p are found from equa- 

tions (A4) and (A5). 

Components of observer-vehicle thrust forces resolved in the target local vertical 
(that is, Tx, Ty, and Tg) are defined as 


T x = a ll T X,b + a 12 T Y,b + a 13 T Z,b 

(All) 

Ty = a 2 iT x>b + a 2 2T Y)b + a 23 T z>b 

(A12) 

T Z = a 31 T X,b + a 32 T Y,b + a 33 T Z,b 

(A13) 

where T X)b , Ty >b , and Tg b are body-axis force summations and ay 
cosines relating the observer vehicle body axes to the target local vertical, 
force summation, for the Gemini thruster configuration, can be written as 

are direction 
The body- 

T X,b = T X+ K 1 + t X- k 2 + ( t Y+ + t Y- + T Z+ + t Z-) k 3 

(A14) 

T Y,b = ( t Y+ " t Y-) k 4 + T x]y- “ 

(A15) 

T z,b = ( T z+ - t z-) k 2 + ( t y- ■ t y+) k 5 + T e+ “ T e- 

(A16) 
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where T x± , T Y± , T z± , T^ ± , and T are "on-off" signals of appropriate magni- 
tude, and the constants Kj are defined as 

Ki = cos 7.6° 

K 2 = cos 26° 

K 3 = sin 26° 

K 4 = cos 26° cos 5° 

K 5 = cos 26° sin 5° 

No provision has been made for simulating individual jet failures. The direction cosines 
a^j are computed as follows: 


a ll = ^11 s * n + ^13 cos 0 s 

(A17) 

a 12 = b 21 sin e s + b 23 cos e s 

(A18) 

a 13 = b 31 sin 6 s + b 33 cos 6 s 

(A19) 

a 31 = -bn cos 0 S + bjg sin Q s 

(A20) 

a 32 = _b 21 cos e s + b 23 sin e s 

(A21) 

a 33 = -b 31 cos e s + b 33 sin e s 

(A22) 

a 21 = b 12 

(A23) 

a 22 = b 22 

(A24) 

a 23 = b 32 

(A25) 


where 0 g has been defined previously (fig. 11 ) and b^j are direction cosines relating 
the observer vehicle body axes to inertial space. The bjj elements, specified as a 
function of four parameters (quaternion elements), are written as 
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bn = 2(a 2 + d 2 )- 1 

(A26) 

bj 2 = 2(ab - cd) 

(A27) 

bi 3 = 2(ac + db) 

(A28) 

b 2 i = 2(ab + cd) 

(A29) 

b22 = 2 (b 2 + d 2 ) - 1 

(A30) 

t»23 = 2(cb - ad) 

(A31) 


b 3 i = 2(ac - bd) 

(A3 2) 

bs 2 = 2(cb + ad) 

(A3 3) 

b 33 = 2 ( c2 + d 2 ) - 1 

(A34) 

The parameters a, b, c, and d are bounded continuous functions of time and can be 
generated by knowing the components of the observer vehicle inertial angular velocity 
resolved in the body axes. The x, y, and z angular velocity components are defined 
as p, q, and r, respectively. Expressed in terms of p, q, and r, the differential 
equations defining the rotation parameters are 

a = ^(-pd - qc + rb) - K £ ae 

(A35) 

b = i^pc - qd - ra) - Kg be 

(A3 6) 

c = i(-pb + qa - rd) - Kg ce 
2 ' 

(A3 7) 

d = + qb + rc) - K e de 

(A3 8) 

where e is defined as 


e = a 2 + b 2 + c 2 + d 2 - 1 

(A3 9) 


and Kg is a gain constant which is determined empirically on the computer. 
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The angular momentum vector is related to the angular velocity by a linear trans- 
formation. By using this basic definition, the body angular rates are computed as 
follows: 

p = Iii'Lj + I 12 L 2 + 113'Lg 

(A40). 

q = l2l' L l + l22* L 2 + *23 l 3 

(A41) 

r = I 31 T L 1 + 

(A42) 

where Ijj' are the elements of the inverse inertia matrix, 
(symmetric matrix). 

and by definition, Iy' = 1^' 

The components of angular momentum Li, L 2 , and Lg can be generated as a 
function of time by knowing the external moments acting on the observer vehicle. The 
differential equations defining the angular momentum components are 

Li = rL 2 - qLg + Mi 

(A43) 

L 2 = pLg - rLi + M 2 

(A44) 

L 3 = qLj - pLg + Mg 

(A45) 


where Mj, M 2 , and Mg are external moment summations resulting from attitude 
jets located on the observer vehicle. The moments are computed from the following 
equations 


Mj = (: 

m i)t^ + - ( m i)t^_ + 

( m i)t 0 . - ( m i)t 0+ 

+ 

NtV - ( m i)t^._ 

+ 

( M i)t x+ + ( M i)t x . 

+ ( m i)t y+ - ( m i)t y . 

+ 

Mt z+ - ( m i)t z . 

m 2 = ( 

m 2)t^ + - ( M 2)t 0 . + 

( m 2)t x _ “ ( M 2 )t x+ 

- 

( M 2)T y+ + ( M ?)T y J 

+ 

Wt z . - ( M 2 )t z+ 




M 3 = .( 

m s)t^ + - Nt^. + 

( m 3)t x _ “ ( M 3 )t x+ 

1 

+ ; 

( M o)t y+ . - ( M 3)T y _ 

- 

— 

( M 3)t z+ + ( M 3 )t z _ 





(A46) 


(A47) 


(A48) 


Standard convention has been adopted for positive yaw, pitch, and roll. 
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The Agena target vehicle is assumed to be passive; however, tumbling capability 
in three degrees of freedom was simulated. The equations defining the target rotational 
motion are identical in form to those previously given for the observer vehicle, with the 
exception = M 2 = M 3 = 0. The necessary equations are summarized below: 

The transformations relating target body axes to inertial space are 


b n ' = 2(e 2 + f 2 ) - 1 (A49) 

t>12' = 2 (eg - hf ) ( A5 °) 

b 13 ’ = 2(eh + fg) (A51) 

b 21 ’ = 2(eg + hf) (A52) 

b 22 ’ = 2(g 2 + f 2 ) - 1 (A53) 

b 2 3 ' = 2(hg - ef) (A54) 

b 31 ’ = 2(eh - gf) (A55) 

b 32 ’ = 2(hg + ef) (A56) 

b 33 ' = 2<h2 + f2) _ 1 (A57) 

The rotation parameters e, f, g, and h are computed from 

e = |(~Pt f - <lt h + r ts) ' K e ee t (A58) 

g = |(p t h - q t f - r t e) - K e ge t (A59) 

h = |(-p t g + 9t e “ r t f ) " K g he t (A60) 

f = |(Pte + q t g + r t h) - Kgfe t (A61) 
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where = e^ + g2 + h^ + f2 - 1 and K e is constant. The target vehicle body rates 
are given as 


p t = Jll’Hj + J 12 'H 2 + J 13 H 3 

(A62) 

<lt = J 2 l' H l + j 22' H 2 + j 23 h 3 

(A63) 

r t = J 3l' H l + J32' h 2 + j 33’ h 3 

(A64) 

where Jjj* are elements of the inverse inertia matrix, and Hj are components of the 
target vehicle angular momenta. The angular momentum components are computed as 
follows: 

Hi = r t H 2 - q t H 3 

(A65) 

R 2 = PtH 3 - r t H x 

(A 66 ) 

H 3 = QtHi - P t H 2 

(A67) 

Simulation Hardware Equations 

This section is primarily devoted to the geometric relationships existing between 
the target image and the pilot’ s line of sight. The vector represented by the directed 
line segment from the pilot's eye to the center of gravity of the target vehicle defines the 
space line of sight (los). The components of the vector R = -x,-y,-z resolved in the 
observer vehicle body axes are given as 

R X,o = "( a ll x + a 2iy + a 31 z ) 

(A 68 ) 

r Y,o = "( a 12 x + a 22y + a 32 z ) 

(A69) 

R Z,o = "( a 13 x + a 23y + a 33 z ) 

(A70) 

Combining equations (A 68 ) to (A70) with the pilot's eye offsets yields 


Ax = R X,o - x p 

(A71) 

A y = r y,o - y P 

(A72) 

Az = Rz,o “ z p 

(A73) 


66 



APPENDIX A 


The angles a and /3 which specify the direction of the line of sight (los) with 
respect to the observer are defined as for mirror drives as 


a = tan"l ^ 
Ax 


)3 = tan - l 


-Az 


[(Ax) 2 + (Ay) 2 ] 1 / 2 


(A74) 

(A75) 


The relative range ^that is, | Ri os |^ which is used to drive the range bed is com- 
puted as follows: 


R los = + [^ x ) 2 + (Ay) 2 + (Az) 2 


2ll/ 2 


(A76) 


The direction cosines relating the target body axes to the line-of-sight axis 
system are now determined. This transformation is represented by the matrix 


D = 


0rL M y Mz> w hich is obtained by applying three successive simple rotations 


about the specified axes. Alternately, the matrix D can be expressed in terms of 
known information; that is, 


(A77) 


»- 0 y H z bb' t 

where [ a and My are s i m Pl e rotations about the Z-observer body axis and the new 

Y-axis, respectively, and the matrices B and B ,T are defined by the matrix elements 
given by equations (A26) to (A34) and equations (A49) to (A57). The superscript T 
denotes the matrix transpose operation. Expanding equation (A77) yields 


l 31 


^21 = “111 s i n a + ^21 cos a 
d22 = -fi2 s i n a + ^22 cos a 
d 23 = -fj 3 sin a + f 23 cos a 


sin (3 

(A78) 

sin j3 

(A79) 

sin /3 

(A80) 


(A81) 


(A82) 


(A83) 
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dsi = fjj cos a sin 0 + f 2 j sin a sin 0 + cos 0 (A84) 

^32 = *12 cos a sbl P + * 22 s * n a s * n P + * 32 cos P (A85) 

d 33 = f ^3 cos a sin 0 + f 23 sin a sin 0 + cos 0 (A86) 

where fjj coefficients are elements of the matrix BB'^\ The matrix elements of B 
and B' are known functions of time, and hence the fjj elements are uniquely specified. 
These elements are computed from the following equations: 


f ll = b ll b ll' + b 12 b 12' + b l 3 b l3 ' 

(A87) 

f 12 = b ll b 2l' + b 12 b 22* + b 13 b 23' 

(A88) 

f 1 3 = b ll b 3l’ + b I2 b 32’ + b l 3 b 33 T 

(A89) 

f 21 = b 21 b ll’ + b 22 b 12' + b 23 b 13* 

(A90) 

f 22 = b 21 b 2l’ + b 22 b 22’ + b 23 b 23' 

(A91) 

f 23 = b 21 b 3l' + b 22 b 32' + b 23 b 33’ 

(A92) 

f 31 = b 31 b ll f + b 32 b 12’ + b 33 b l3’ 

(A93) 

f 32 = b 31 b 2l' + b 32 b 22* + b 33 b 23' 

(A94) 

f 33 = b 31 b 3l' + b 32 b 32* + b 33 b 33' 

(A95) 


Cockpit Instrument and Auxiliary Equations 

Target azimuth and elevation with respect to the local vertical axis system is 
given by 


X = tan -1 ^ (A96) 

6 = tan -1 -2- (A97) 

K xy 
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where 


R xy = -(x cos X + y sin X) (A98) 

Range rate is computed for display purposes and is given by 


R = 


xx + yy + zz 
Kos| 


(A99) 


where |Rios| is defined by equation (A7 6). The components of Av are computed from 
the following definitions: 


AV X 


AVy 


A V z 



T X, b 

T Y, b 

T Z, b 


dT 

dr 

dT 


(A100) 

(A101) 

(A102) 


In these equations, the total mass m is assumed to be constant. The mass changes due 
to translational and attitude control jets, respectively, are given by 


4mtrans = 

Thus, the percent fuel used 
puted as follows: 


f 

J 0 


T X+ + | T X- + T Y+ + T Y _ + T Z+ + T Z J dr (A103) 




+ T 0+ + T e _ + T^ + + T 0 _|)dr (A104) 


by the translational and attitude control systems are com- 


F trans _ 


^ m trans 


100 


F att 


100 

F o 


where F Q is total fuel initially in kilograms. 


(A105) 

(A106) 
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DERIVATION OF GREEN'S FUNCTION TECHNIQUE 


This appendix is intended to give a comprehensive derivation of the Green's func- 
tion technique used in this paper. The motivation for this derivation is that the equation 



(Bl) 


for the values of £ and w s for this problem and related problems is difficult to 
solve by using standard numerical techniques. The difficulty in solving this equation is 
not the forcing function p(r), but is found to be the differential operator on the left-hand 
side of equation (Bl). The approach used in this paper was to find a different linear 
operator, generically different from the differential operator, with the hope that a more 
adequate numerical approximation exists for this new operator. Such an operator was 
found which cast equation (Bl) as an integral equation involving Green's function. 

The Green's function of equation (Bl) satisfies the equation (ref. 19) 

G(t - r) = w g 2s(t - t) (B2) 

The Green's function also satisfies boundary conditions. These values at the boundaries 
must be chosen so that a solution to equation (B2) exists and the integral representation 
of equation (Bl) has a simple form. 

Multiplying equation (Bl) by G(t - r) and equation (B2) by Pg( T ) and subtracting 
the two results gives 



d 9 
W S dt + 


- T >^ + 2 ?g“g IF + w g 2 )%< T > - + 2 h“s 2 Tt 


+ w g 2 ]G(t - r) 


= w g 2 |G(t - t) p(r) - Pg(r) 6 (t - rjj (B3) 
Integrating equation (B3) on the interval [ 5 , t + ej results in 
w g 2 j* |G(t - t) P (t) - Pg(r) 6 (t - t)] dT 


-£ 


0 
it+e 

0 


° (t " T) (^ + + W g 2 ) P g (T) " P g (T) fe + 2? g C °g * + ^g^) 0 ^ “ T > 


dr (B4) 
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Performing the integration of equation (B4) with the delta function yields 
p g (.) - Oft - r) P(T) dr ♦ y o t+£ <jp g(T )gG(t - ri] - G(t - r) 

2 ?e P t+€ 




>dr 


o>, 


g 


C v g(T) fe G(t ' T ] ' G(t ' T) &% (T ]r T 


(B5) 


Consider the last two integrals of equation (B5) . Note that 


Then 


!«' - » - - T > 


(B6) 


§* [ p g (T) S G(t ■ T i| ■ G(t ■ T) [if p g (T) ]| dT 

= ■ C { p « (T) [s G<t ' T >] + G(t ' T) [l? p g <T | 

= ■ C ^l? g(T) g(i ■ dT 


dT 


= -Pg(T) G(t - r) 


t+€ 

0 


(B7) 


The second integral of equation (B5) is evaluated by integrating by parts. Consider 


r*t+€ 

\ G(t - T) 
J 0 


d d T 2 p e <T)l 


dr = G(t - r) -£-p(r) 


* +e pt+e H H 

0 -} 0 (») 


Integrating again by parts yields 
^t+€ 


ni+t 

\ G(t - T) 
J 0 




t+€ p- 


dT = C (t - T) g- Pg (T ) 

+ C p * (T) & G(t ' T i 


Tr G<‘ - T) 


Pg(T) 


t+e 

0 


dT 


(B9) 


By substitution of equation (B6) and the identity 
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i|G(t-r)=-^_G(t-r) 


dt 2 


dr^ 


equation (B9) becomes 


nt+€ 

G(t - r) 

J 0 


Pp-( T ) 


dr2 » 


d T = G( t - T )[|-p g (r]| o + [|G(t-r^ 

+ r p ^ <T) iS G(t ■ r i 


p g ( T ) 


t+e 

0 


Substituting these results (eqs. (B7) and (Bll)) into equation (B5) yields 
->t +6 

G(t - t) p(r) dr - - 

"'g 


% (t) = c G(t ■ t) p(r) dT ■ ^ r ■ t) ^ Pg(r) + [^° (t ■ pg(T t 


|t+€ 


- 77 S Pg( T ) G (t - t) 

U) g 6 


|t+€ 

0 


(BIO) 


(Bll) 


(B12) 


Evaluating the algebraic expressions at their limits and grouping yields 


n t+e 

P g (t) = G(t - t) p(r) dT 

G(-e) p g (t + e) + G(-e) p g (t + e) + 2? g G(-e) p g (t + e)^ 

G(t) Pg(0) + G(t) p g (0)] + 2? g G(t) p g (o| 


CUgW 


w g K 


(B13) 


with the convention that f(x) = lim Therefore a solution to equation (B2) is needed 

y-+x Q y 

which satisfies the following conditions: 


■y 


G(t) = 0 


(t ^0) 


lim G(-t) = 0 j 

t^O l 


6(0) = co g 2 

lim G(t) = lim G(t) = 0 

t->oo 


(B14) 
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Taking the limit of equation (B13) as e goes to zero yields 


p g (t) = y G(t - t) p(r) dT + 


J_G(t) + ^G(t) 


COc 


'g~ s 

Adding the further restriction on the Green's function yields 


p g (0) + ^2 G(t) % (0) (B15) 


G(0) + 2? g w g G(0) = 0 (B16) 

The boundary conditions given by equations (B14) and (B16) are needed for the solution to 
have a simple form (eq. (B15)) which satisfies the initial conditions of Pg(t). A solution 
to equation (B2) must be found which satisfies these boundary conditions, that is, if such 
a solution even exists. 

The existence of the Green T s function is shown by finding a solution to equation (B2) 
which satisifes equation (B14) and equation (B16). Consider the equation 


(S +2 5i“g$ + “ g 2 ) G <« = “ g 2 


Taking the Laplace transform s of equation (B17) (ref. 19) yields 


(B17) 


s 2 + 2£ g co g s + co g 2 ^G*(s) = 


CO 


g 


(B18) 


from which G*(s) is found, as expected, to be the transfer function of equation (Bl) 

2 


G*(s) = 


co„ 


s2 + 2? g co g s + co g 2 


(B19) 


The inverse Laplace of G*(s) can be found in standard tables to be 


G(t) = — £-e - ^ sin rt 
2\l/2 


(t g 0) (B20) 


/ 2 \ 1/2 

where £ = ?gCo g and T = C0g(l ?g J . The first two derivatives of equation (B20) 

G(t) = -|G(t) + co g 2 e“^ cos rt (B21) 


are 


G(t) = -|G(t) - ijcOg 2 e & cos rt - T 2 G(t) 


(B22) 
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for t £ 0 and by equations (B14) 


G(t) = G(t) = 0 (t < 0) 

Hence, the following boundary conditions hold: 


G(0) = 0 


lim G(-t) = 0 
t^O 

G(0) = w g 2 

lim G(t) = lim G(t) = 0 

t->00 t-Z'OO 

From equation (B22) 

G(0) = -|G(0) - ^ g 2 = -2£ g w g G(0) 

Thus, the Green's function of equation (B20) meets all the requirements of equations (B14) 
and (B16). 

In determining equation (B15), no assumptions were made of the form of p(r). 
Equation (B15) is a valid solution even when p(r) is a nonlinear function of p g . 

The remainder of this appendix is centered on the evaluation of the integral. Nor- 
mally, initial conditions are zero which is consistent with the usage in the text. Initial 
conditions can be incorporated in the resulting integration scheme; this procedure is 
shown later in the appendix. Equation (B15) with zero initial conditions is 

p g (t) = J G(t - t) p(r) dr (B23) 

which is not in a useful form. The solution would require storage of the forcing function 
from the beginning of the problem and the evaluation of a rather complex integral over 
increasingly large ranges as the problem progressed. For an efficient integration of 
equation (B23), a form needed to be found so that only local information of the forcing 
function and values of the integral are needed to propagate the solution. To find such a 
form, consider p (t) to be known at time t and p(r) specified on the interval (t, t + h) 

o 

and look at equation (B23) at time t + h; then, 


/-.t+h 

p g (t + h) = J G(t + h - t) p(r) dr (B24) 
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Note that 


G(t + h - r) = -r_[ G t(h) G(t - r) + G(h) G^t - t)| 


(B25) 


where 


g 


G^(t) = — £-e~^ cos rt 


-«t+h 


Equation (B24) is then 

-p 4 . pt+h .p pt+n + 

p„(t + h) = — - G'(h) l G(t - t) p(r) dr + — ~ G(h) \ G'(t - r) p(r) dr 
& J 0 w g 2 ^0 

pt pt+h 

\ G(t - t) p(r) dr + \ G(t - r) p(r) dT 
Jq J t 

(* G^(t - t) p(r) dT + f G^(t - t) p(r) dr 
J 0 J t 


= J^G t (h) 


+ G(h) 

V 


(B26) 


where equation (B26) has the form of an integration over the last interval with the addi- 
tion of previously computed information provided the auxiliary computation 


A(t) = C* G f (t - r) p(r) dr 
Jo 


(B27) 


is made. 


The computation of equation (B27) can be placed in a similar form by using the 
same process. Consider 

-.t+h 


pL+n . 

A(t + h) = \ G'(t + h - r) p(r) dr 
J 0 


where 


G t (t + h - r) = ^G t (h) G t (t - r) - -^G(h) G(t - r) 

w g 


u » g 2 


(B28) 


(B29) 


Equation (B28) can then be written 
A(t +h) = T G t( h) 


u ) g 2 


G(h) 


f G^(t - t) p(r) dr + C G^(t - r) p(rl dT 

Jo Jt 

pt pt+h 

\ G(t - t) p(t) dT + \ G(t - t) p(t) dT 

Jq Jt 


(B30) 
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Substituting equations (B27) and (B23) into equations (B26) and (B30) and writing the 
result as a matrix product results in 


— — 
p g (t + h) 

r 

G f (h) 

G(h) 

pt+h 

p g (t) + J G(t - t) p(r) dr 

A(t + h) 

co 2 
g 

-G(h) 

G f (h) 

^t+h + 

A(t) + J G' (t - t) p(r) dr 


(B31) 


To this point no approximations or assumptions have been made. The solution is 
valid and almost in a form to be realized on a digital computer. At this point many vari- 
ations of the evaluation of the integrals leading to various numerical techniques can be 
made. An example is the approximation for p(r) on I(t, t + h) of 

p(r) * p(r) 5(t - t) (B32) 


which leads to the primitive recursive filter. (See ref. 5.) In this application p(t) is 
known to the accuracy of the integration scheme in the moment equations where p(r) is 
computed. Hence, p(r) is known to be some polynomial on I(t, t + h): 

p(r) = a + br + cr 2 + . . . fr n (B33) 


Within this framework, equations (B31) can be evaluated exactly (at least within computer 
accuracy). To do this, first change the integration variable by the transformation 


Thus, 


x = t - t (B34) 

dx = dr (B35) 


and equations (B31) become 


p g (t + h) 

r 

G^h) 

G(h) 

P g (t) + j 

->h 

G(-x) p(t + x) dx 

'0 

A(t + h) 

u»g 2 

-G(h) 

G T (h) 

A( t ) + 

h t 

G'(-x) p(t + x) dx 

) 


and equation (B33) becomes 


n 

p(t + x) = ^ ajX^ 

3=0 


(B36) 


(B37) 
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Substituting equation (B37) into equations (B36) yields 






n 

P g (t + h) 

r 

“cog 2 

G t (h) 

G(h) 

p g (t) + 1 a j x i 
j=0 

n 

v — ' 

A(t + h) 

-G(h) 

G^h) 

A(t) + ^ a jlj ' 

_ 




)=° 


where 

r h 

L = \ G(-x)x^ dx 

J Jo 

= C G^(-x)x^ dx 
J Jo 

for j = 0, 1, 2 . . . . Rewriting equation (B21) in terms of G(t) and G^(t) 


G(t) = -£G(t) + rc t (t) 


so that 

G(t) + 2£ g co g G(t) = rG T (t) + |G(t) 


Neglecting the integral part of equation (B15) results in 


p g (t) = 


-^G^t) +-i^G(t) 


Pg (0) + — “2 G(t) p g (0) 


Advancing equation (B41) in time results in 


p„(t + h) = — ^rP 0.(0) G*(t + h) + 


^2*8' 
g 


9 P g (0) + -^p (0) 

CO B COg^ s 


G(t + h) 


= ^p g (0)[pt(h) G*(t) - G(h) G(t[| 


— ^TPg(O) + — ^5-Po.(0) 


CO. 


g 


CO, 


2 g 


G^h) G(t) + G(h) G^t) 


(B38) 


(B39) 


(B40) 


(B41) 


(B42) 
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Grouping terms in equation (B42) gives 


P e (t + h) = -^G^hX 
& oo * 


g 


^Gt(t) + -^G(t) 

w g 


L w g 2 


Pp-(O) + — ^ Pg(0) G(t)' 


oo 


g 


^GChJF-^GXt) - -^G(t) 

OOg * | Wg* 4 00 „ 


g 


p g (0) + r2 p g (0) ° t(t) > 


00 


g 


Comparison of equation (B43) with (B38) yields for p(t) = 0 


A(t) = 


— G*(t) - ~^»G(t) 
V oo g l 


Pg(°) + — ^ Pg(0) G*(t) 

Cc) nr 


The initial conditions to start the integration scheme (eqs. (B38)) are then 


P g (0) = P g (0) 


A(0) = |;Pg(0) + |;Pg(0) 


(B43) 


(B44) 
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TABLE I.- ROOT LOCATIONS FOR EULER PARAMETERS 
AND DIRECTION COSINES 


Parameter 


Euler 

parameters 


Direction 

cosines 


Laplace 

Adams 

Bashforth 

Euler 

Time 

constant 

Frequency 

Time 

constant 

Frequency 

Time 

constant 

Frequency 

oo 

o>/2 

From equation (69) 

From equation (64) 



1 256 
4 h 3 w 4 

2 96 

8 

ho> 2 

o> h 2 w 3 
2 24 

oo 

c 0 

From equation (59) 

From equation (41) 



1 16 
4h3 w 4 

u> + -j^h3<x>3 

2 

ho> 2 

to - 

3 
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TABLE H.- ROUNDOFF AND TRUNCATION ERRORS 


IN TRANSLATIONAL VARIABLES 


Source j 

p 

— 

8 , 

radians 

x, 

meters 

z, 

meters 

Time, 1008 sec 

Double precision Runge-Kutta (4th order) 

Partial double precision Adams-Bashforth (2d order) 
Single precision Adams-Bashforth (2d order) 

-0.00299817 

-.00299817 

-.00299756 

1.1991391 

1.1991391 

1.1990303 

-1179.7356 

-1179.7333 

-1179.5947 

-1381.6586 ! 

-1381.6559 

-1381.3930 


Time, 2004 sec 


Double precision Runge-Kutta (4th order) 

-0.01266422 

2.3510946 

1990.8761 

-1055.6684 

Partial double precision Adams-Bashforth (2d order) 

-.01266422 

2.3510947 

1990.8720 

-1055.6661 

; Single precision Adams-Bashforth (2d order) 

-.01266017 

2.3506772 

1990.1228 

-1055.2903 




Time, 3000 sec 


Double precision Runge-Kutta (4th order) 

-0.01626021 

i 1 

3.4786051 

-2717.4585 

486.8879 

Partial double precision Adams-Bashforth (2d order) 

.01626021 

3.4786052 

-2717.4524 

486.8873 

Single precision Adams-Bashforth (2d order) 

.01625222 

| 3.4777238 

-2715.9020 

486.6915 


Time, 4008 sec 


Double precision Runge-Kutta (4th order) 

0.00346380 

4.6313472 

-148.1033 

1474.5843 

Partial double precision Adams-Bashforth (2d order) 

.00346380 

4.6313472 

-148.1018 

1474.5815 

Single precision Adams-Bashforth (2d order) 

.00346282 

4.6297367 

-147.7157 

1473.6415 


Time, 5004 sec 


Double precision Runge-Kutta (4th order) 

-0.01052772 

5.8052815 

-2703.0378 

681.0652 

Partial double precision Adams-Bashforth (2d order) 

-.01052772 

5.8052815 

-2703.0333 

681.0635 

i Single precision Adams-Bashforth (2d order) 

-.01051593 

5.8027082 

-2701.1182 

680.5329 






Figure 3.- Predicted and computed errors in the norm of the direction cosines for p = q = 

and h = ■— second using Euler integration. 


radian/second 


<D 

U 

u 

<D 

a 


u 

o 

u 

u 

w 


1.0 



Figure 4.- Predicted, computed, and starting errors in the norm of the direction cosines for p = q = r = radian /second 
and h = — second using Adams-Bashforth second-order integration with Euler integration as starting formula. 
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Figure 5.- Predicted and computed errors in norm of the Euler parameters for p = q = 

and h = ~ second using Euler integration. 


: -j= radian/second 


0 . 125 



Figure 6.- Predicted, computed, and starting errors in norm of the Euler parameters for p = q = 
1 


ft 


radian/second 


and h = V 7 second using Adams-Bashforth second-order integration with Euler integration as starting formula. 
16 
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0. 0 o. 125 0. 25 o. 375 

t, sec 

Figure 7.- Gyroscope response (solid) to a step input (dashed) for oig = 125 radians/second and £g = 0.8 with h = 

using the Green's function technique. 



t, sec 

Figure 8.- Gyroscope response (solid) to a ramp input (dashed) for oig = 125 radians/second and ^g = 0.8 with h 

using the Green's function technique. 




0 40 80 120 

t, sec 


Figure 11.- Predicted and computed errors in total angular momentum for Ijj = 500.0 kg-m 2 , |^2 = 878.9 kg-m 2 , 
I 33 = 881.9 kg-m 2 , p = q = r = -]= radian/second and h = ~ second using Euler integration. 



t, sec 


Figure 12.- Predicted, computed, and starting errors in total angular momentum for 1^ = 500.0 kg-m 2 , I 22 = 878.9 kg-m 2 , 
I 33 = 881.9 kg-m 2 , p = q = r = j= radian/second and h = ^ second using Adams-Bashforth second-order integration 
with Euler integration as starting formula. 
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(a) h = 31.25 msec. 


h - 1 5 . 625 msec 



10. OQ 20.00 

t, sec 


(b) h = 15.625 msec. 

Figure 13.- Time-history plot of p for attitude control step input. 


radians/sec radians/sec 



t, sec 

(a) h = 31.25 msec. 



t, sec 


(b) h = 15.625 msec. 

Figure 14.- Time-history plot of q for attitude control step input. 






0.0 0.0625 0.125 0.1875 0.25 

h, sec 

Figure 16.- Global error as function of step size h in trajectory variables at t = 5004 seconds for Euler (E) and Adams-Bashforth ( AB) 
second-order integration over the range ^ second = h =| second on a 27-bit fractional part machine. 
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Figure 18.- Equation flow diagram. 






















National Aeronautics and Space Administration 
Washington, D. C. 20546 


POSTAGE AND FEES 
NATIONAL AERONAUTI 
SPACE ADMINISTRA1 


OFFICIAL BUSINESS 


FIRST CLASS MAIL 



i L I 


. I 1 


POSTMASTER: 


If Undeliverable (Sect 
Postal Manual ) Do N 


The aeronautical and space activities of the United States shall be 
conducted so as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. ’The Administration 
shall provide for the widest practicable and appropriate dissemination 
of inf or ??i at ion concerning its activities and the results thereof.” 

— National Aeronautics and Space Act of 1958 ,, 

/ 

NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technplogy Utilization Reports and Notes, 
and Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION DIVISION 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 



